Identifying the Effect of Public Holidays on Daily Demand for Gas

Summary

To reduce operational costs and to ensure security of supply, gas distribution networks require accurate forecasts of the demand for gas. Among domestic and commercial customers, demand relates primarily to the weather and patterns of life and work. Public holidays have a pronounced effect which often spreads into neighbouring days. We call this spread the ‘proximity effect’. Traditionally, the days over which the proximity effect is felt are prespecified in fixed windows around each holiday, allowing no uncertainty in their identification. We are motivated by an application to modelling daily gas demand in two large British regions. We introduce a novel model which does not fix the days on which the proximity effect is felt. Our approach uses a four-state, non-homogeneous hidden Markov model, with cyclic dynamics, where the classification of days as public holidays is observed, but the assignment of days as ‘pre-holiday’, ‘post-holiday’ or ‘normal’ days is unknown. The number of days to the preceding and succeeding holidays guide transitions between states. We apply Bayesian inference and illustrate the benefit of our modelling approach. A version of the model is now being used by one of the UK's regional distribution networks.

Calendar effects, Forecasting, Gas consumption, Hidden Markov model, Time series

1 Introduction

The energy sector in the UK is changing. To comply with the UK Climate Change Act 2008, greenhouse gas emissions must be reduced to 80% of their 1990 levels by the year 2050 (Her Majesty's Parliament, 2008). Pragmatically, these targets must be met in a manner which does not compromise the provision of affordable prices to consumers or the competitiveness of UK industry. At present, natural gas provides a comparatively low cost source of energy which, unlike some renewable energy sources such as wind power, can also offer reliability in supply and storability (Chu and Majumdar, 2012). In the future, decarbonization of the gas network is likely to involve connecting carbon neutral biomethane plants to the transmission system or converting the distribution network to transport hydrogen rather than natural gas (Dodds and McDowall, 2013). Therefore, both now and over the years to come, gas has a vital role to play in the energy mix and it is more important than ever that the gas distribution network operates as efficiently as possible.

National Grid is the sole owner and operator of the gas transmission infrastructure in the UK. Gas in the national transmission system leaves the network at high pressure at 49 points across the country. After being odorized for safety it is transported, ultimately at lower pressure, through eight regional distribution networks to individual customers. National Grid works closely with the regional distribution networks to ensure that the local supply of gas meets the demand at all times.

Demand forecasts, over a range of horizons, are required by both National Grid and the distribution companies for reasons including safety and security of supply and investment and operational planning (National Grid, 2016). The method that is used by distribution companies involves forecasts ${\widehat{A}}_{u,j}$ for the annual demand A_{u, j} in year u from a subset j of customers, j = 1,…, J. These are supported by a report, produced by an external company, taking into account economic and other external factors. Forecasts for demand ${\tilde{Y}}_{t,j}$ for subset j on day t are then rescaled to match this annual forecast; see National Grid (2016). Thus what is required is inference about a set of scale factors exp (k_{t, j}) for the days of the year so that the forecast mean of ${\tilde{Y}}_{t,j}$⁠, for example, will be $E\{\exp (k_{t,j})\}{\widehat{A}}_{u,j}$⁠, since the scale factors may reasonably be treated as independent of the annual demand level.

In this paper we describe a model for demand using the transformation $Y_{t,j}=\ln ({\tilde{Y}}_{t,j})=c_{u(t),j}+k_{t,j}+e_{t,j}$ where u(t) is the year in which day t falls. Here e_{t, j} represents random fluctuation and the value of the additive constant c_{u(t), j} can be adjusted in forecasts to match ${\widehat{A}}_{u(t),j}$⁠, without affecting the log-scale-factors k_{t, j}. It is these factors, incorporating weather, seasonal and calendar effects, about which inference is required.

There is a large body of work in the scientific literature concerned with modelling and forecasting the demand for natural gas; see, for example, Soldo (2012) for an extensive review. In broad terms, consumers of gas can be divided into three types: residential, commercial and industrial. We focus on the first two of these groups in this paper. Although economic factors like gas prices and national income are important drivers of gas consumption by big industrial users, their effects are generally less important for residential and commercial customers for whom gas is predominantly used for heating (including water) and cooking. As a result, the demand by these groups is strongly related to the weather and patterns of life and work. Models for residential and commercial gas consumption therefore generally allow for weather-related predictors and seasonal and calendar effects. Sometimes models also incorporate interactions between the two to allow the effect of the weather and, in particular, temperature to vary with periodic changes in fixed or timed heating schedules. After modelling the weather, seasonal and calendar effects, it is common for remaining errors to exhibit considerable auto-correlation. This is generally modelled directly, for instance by assuming auto-regressive moving average models for the residuals (Lyness, 1984; Akouemo and Povinelli, 2016).

Weather-related predictors in models for gas consumption generally include primitive variables, like temperature and wind speed, or derived variables, constructed to have a simple relationship with demand. In the literature, considerable attention has been devoted to modelling the non-linear effect of temperature which manifests at low and high temperatures when, for example, many customers switch their heating on or off. In some cases this is achieved by fitting non-linear regression models (for example see Brabec et al. (2009) and Gascón and Sánchez-Úbeda (2018)) whereas, in others, the effect is captured through bespoke, derived variables such as ‘heating degree days’ (e.g. Akouemo and Povinelli (2016)). The latter measure the temperature difference from some fixed threshold but adopt 0-values when the temperature exceeds that threshold, so that the temperature has no effect when heating is not needed.

Seasonal and calendar effects often include a smoothed effect for the day of the year (e.g. Brabec et al. (2009) and Gascón and Sánchez-Úbeda (2018)) and fixed effects to represent the day of the week. These are needed because gas demand follows a weekly cycle, with clear differences between weekdays and weekends (Lyness, 1984). Similarly, because demand is affected markedly by public holidays, models usually include a holiday factor (Brabec et al., 2015) or treat public holidays like days of the weekend (e.g. Charlton and Singleton (2014)). In addition to effects on the days of the holidays themselves, models for energy consumption sometimes include a protracted effect which extends into neighbouring days. This allows for the influence of changes in cooking and travel arrangements, and for the commercial and industrial slow-down that typically occurs around holidays. For example, in their model for gas consumption in the Czech Republic, Brabec et al. (2009) included Christmas and Easter effects which apply over a fixed range of consecutive days around Christmas Day and Good Friday. Taking a different approach, Brabec et al. (2010) defined a ‘daytype’ factor with five levels corresponding to different classifications of the current, previous and following days as working days or otherwise. This allows public holidays and neighbouring days to have different effects depending on where they fall in the week. In their models for electricity consumption, Hor et al. (2006) and Pardo et al. (2002) allowed demand to differ on proximity days within a short fixed window around each holiday.

Some public holidays, particularly those around Christmas and often Easter, occur during periods where the demand for gas is at its highest. At such times, it is essential to have accurate forecasts and a correct quantification of uncertainty. Yet, in terms of allowing for a proximity effect, an aspect that is common to all of the approaches above is that the dates of the proximity days are fixed, often arbitrarily.

In this paper, we are motivated by an application involving daily gas consumption in the two large geographical areas in northern England served by the regional operator Northern Gas Networks (NGN). We propose a novel model for daily gas demand which, to our knowledge, is the first approach which does not fix the days on which the proximity effect of a public holiday is felt. The results of a preliminary version of this model are already being used by the company in its annual medium-term forecasting exercise. Our approach is based on a four-state, non-homogeneous hidden Markov model with cyclic dynamics. In this model the classification of days as public holidays is observed, but the assignment of days as ‘pre-holiday’, ‘post-holiday’ or ‘normal’ days is unknown but guided by the number of days from the preceding and to succeeding public holidays. We allow for auto- and cross-correlation in the bivariate time series by modelling the logarithm of gas demand in the two regions, conditionally on the states, using a bivariate auto-regression of order 1, with a symmetric auto-regressive coefficient matrix, that allows stationarity to be imposed through simple constraints on the parameter space. Taking a Bayesian approach to inference, we use a hierarchical structure to encapsulate structural prior information about similarities between the two regions. In our application, we illustrate the benefit of allowing for a proximity effect of unknown magnitude and duration.

The remainder of this paper is structured as follows. The data for our analysis are described in Section 2. In Section 3 we propose our model and describe it in detail. Section 4 gives the structure of our prior distribution and Section 5 discusses computation of the posterior distribution. We give the results of applying our model and inferential procedures in Section 6 and draw some conclusions in Section 7.

2 The data

NGN, one of the eight regional distribution networks in the UK, is responsible for gas distribution to 2.7 million homes and businesses across a 25000-km² region in the north-east of England, northern Cumbria and Yorkshire.

In the UK, medium- and long-term demand forecasts are produced for each of 13 local distribution zones (LDZs). NGN is responsible for generating the demand forecasts for the neighbouring Northern, NO, and North-East, NE, LDZs. The former covers northern Cumbria and the north-east of England and the latter encompasses much of Yorkshire. Within each LDZ, every point at which gas is taken from the network by supply pipe to a consumer is categorized according to its load band, which is based on how much gas it uses. Depending on whether or not meters are read daily, load bands are classified as daily metered or non-daily metered (NDM). Daily metered load bands comprise large industrial premises which typically have the highest demand for gas. NDM load bands are subdivided into four categories, roughly containing domestic, commercial and small and medium-sized industrial customers. In this paper we focus on load bands 1–3, the definitions of which are shown in Table 1.

Table 1

Classification of NDM load bands 1–3

Index	Load band	Example
	(MWh year−1)
1	0–73	A single house
2	73–732	A large block of flats or commercial premises
3	732–5860	Small industrial premises

For both LDZs and each of the three NDM load bands, daily gas consumption data are available for the 9-year period from January 1st, 2008, to February 18th, 2017. The accompanying weather data take the form of a derived variable called the composite weather variable (CWV), defined by National Grid, which takes a single daily value for each LDZ. The CWV is based primarily on air temperature. However, to strengthen the linear association between this temperature-based variable and gas consumption, its construction accommodates adjustments at high and low temperatures and allowance for the effects of other weather variables, such as wind speed. Further details can be found in National Grid (2016).

3 A non-homogeneous hidden Markov model

Denote by ${\tilde{Y}}_{t,j}$ the gas demand in tenths of a gigawatt-hour pertaining to a given NDM load band on day t in region j, where j = 1 and j = 2 correspond to the Northern NO and North-East NE LDZs respectively. Rather than modelling the raw demand, we instead work on a log-scale, defining $Y_{t,j}=\ln ({\tilde{Y}}_{t,j})$ and Y_t = (Y_t,1, Y_t,2)^′ for t = 1,…, T. This helps to make the variance more stable across seasons and gives fixed effects in our additive model a multiplicative effect on the original scale. The CWV on day t is denoted by w_t = (w_t,1, w_t,2)^′. We introduce further explanatory variables $n_t\in {\mathbb{N}}_0=\{0,1,2,\dots \}$ and $p_t\in {\mathbb{N}}_0$ which indicate the number of days to the next and since the previous public holiday respectively. Clearly day t is a public holiday, known colloquially in England as a bank holiday, if and only if $n_t=p_t=0$⁠. Finally, we introduce the categorical explanatory variable $r_t\in \{1,2,3\}$ which indicates the ‘type’ of the nearest public holiday according to the calendar season in which it falls. The labelling is given in Table 2. Note that if a holiday, such as Christmas Day, falls at the weekend, it is a neighbouring weekday that will be observed as the corresponding bank holiday.

Table 2

Categorization of UK public holidays†

Type	Name	Days
1	‘Easter’	Good Friday and Easter Monday
2	‘Other’	May day, spring bank holiday, summer bank holiday and one-off holidays such as the Queen's Jubilee
3	‘Christmas’	Christmas Day, Boxing Day and New Year's Day

†

If day t is equidistant between two holidays of different types, the type r_t of the later holiday is used.

As discussed in Section 1, public holidays often have a protracted effect on gas consumption, which extends into neighbouring days. The simplest way of modelling this effect would be to classify days within a fixed window around each public holiday deterministically as proximity days. However, this approach is inflexible, relying on an arbitrary decision about the existence, size and position of the windows. Instead, we allow the number of proximity days on either side of each public holiday to be unknown. This is achieved by introducing a discrete-valued stochastic process {S_t:t = 0,1,…, T} with $S_t\in {\mathcal{S}}_s=\{1,2,3,4\}$ in which state 2 can apply to public holidays only and is observable, whereas states 1, 3 and 4 are ‘hidden’ and cannot be observed. Of these, states 1 and 3 allow for pre-holiday and post-holiday proximity days, whereas state 4 is a baseline for normal days. To support this we define a distribution over the states so that state 1 can be entered only from state 4 and left via state 2, whereas state 3 can be entered only from state 2.

The directed acyclic graph in Fig. 1 illustrates the assumptions of conditional independence in our hidden Markov model. A non-homogeneous Markov chain for the states will be described in Section 3.1, whereas the model for gas demand, conditional on the states, will be discussed in Section 3.2.

Fig. 1

Directed acyclic graph illustrating the factorization of the joint density of the observations and the hidden states conditional on the time series of explanatory variables $(n_t,p_t,r_t,{\mathbf{w}}_t)$⁠, for t = 1,2,…: in a directed acyclic graph random variables are represented by nodes; these are connected by directed arrows which indicate the order of conditioning when factorizing their joint probability density

3.1 Modelling the states

The challenge in modelling the joint distribution of the states arises because state 2, representing public holidays, is observable and the dates of public holidays are always known in advance. In the hidden Markov model framework, we can build this information into our prior distribution for the states S_t in two ways. The more direct approach would be for the states to evolve according to a non-homogeneous Markov chain with transition probabilities that vary over time. Transition into state 2 can occur with certainty if day t is a public holiday, and transition into the pre-holiday state (state 1) can become more likely in the days leading to a holiday. Alternatively we can specify the joint prior distribution for the states by updating an ‘initial’ prior, representing a homogeneous Markov chain, by conditioning on the observation that S_t = 2 if day t is a public holiday or S_t ≠ 2 otherwise, for t = 0,…, T. Of course, the resulting process is no longer Markovian.

We choose to adopt the more direct approach, based on a non-homogeneous Markov chain, and model the transition probabilities as functions of the number of days to the next (⁠$n_t$⁠), and since the previous (p_t), public holiday. This approach allows a flexible, non-geometric distribution for the pre- and post-holiday state sojourn times which could, for example, rapidly decay after 2 days, in keeping with views from the literature (see Section 1) and the expert judgement of engineers at NGN (see Section 3.1.2). Moreover, use of covariates $n_t$ and p_t in the non-homogeneous model enables straightforward experimentation with the joint prior induced for the states, which aids in constructing a distribution with a sensible concentration of prior mass (see Section 6.1).

3.1.1 Public holidays

Suppose that the states S_t follow a first-order, non-homogeneous Markov chain with

$$\text{Pr}(S_t=k|S_{t-1}=j,n_t,p_t)={\lambda }_{j,k}(n_t,p_t),(j,k)\in {\mathcal{S}}_s^2,$$

for t = 1,2,…, T, and initial distribution

$$\text{Pr}(S_0=k|n_0,p_0)=l_k(n_0,p_0),k\in {\mathcal{S}}_s.$$

Since state 2 is observable, it follows that, if day 0 is a public holiday, then n₀ = p₀ = 0 and l₂(0,0) = 1. Similarly, for any $j\in {\mathcal{S}}_s$⁠, if day t is a public holiday, then $n_t=p_t=0$ and λ_j,2(0,0) = 1, as depicted in Table 3.

Table 3

Transition matrix for the state process when day t is (a) not a public holiday and (b) a public holiday

St − 1	St
	1	2	3	4
(a) Not public holiday
1	1	0	0	0
2	0	0	λ2, 3(·)	1 − λ2, 3(·)
3	0	0	1 − λ3, 4(·)	λ3, 4(·)
4	λ4, 1(·)	0	0	1 − λ4, 1(·)
(b) Public holiday
1	0	1	0	0
2	0	1	0	0
3	0	1	0	0
4	0	1	0	0

3.1.2 Other days

In the remainder of this section we consider the initial distribution $\text{Pr}(S_0=k|n_0,p_0)$ and transition probabilities $\text{Pr}(S_t=k|S_{t-1}=j,n_t,p_t)$ characterizing the state process if day t is not a public holiday. In such cases $(n_t,p_t)\in {\mathbb{N}}_0^2\backslash \{(0,0)\}$ and state 2 cannot be assigned. Dropping the t-subscript for brevity, it follows that l₂(n, p) = 0 and λ_j,2(n, p) = 0 for all $j\in {\mathcal{S}}_s$⁠. To complete our initial distribution we assign $l_k(n,p)=\frac{1}{3}$ for $k\in {\mathcal{S}}_s\backslash \{2\}$⁠. Thereafter, we impose cyclic dynamics on the evolution of the state by assuming that all the remaining transition probabilities λ_{j, k}(n, p) for n > 0, where $j\in {\mathcal{S}}_s$ and $k\in {\mathcal{S}}_s\backslash \{2\}$⁠, are 0 except the transitions from the normal state to the pre-holiday state, λ_{4, 1}(n, p), transitions from the holiday state to the post-holiday or normal states, λ_{2, 3}(n, p) and λ_{2, 4}(n, p), transitions from the post-holiday state to the normal state, λ_{3, 4}(n, p), and the self-transitions λ_{1, 1}(n, p), λ_{3, 3}(n, p) and λ_{4, 4}(n, p). Necessarily this implies that λ_{2, 4}(n, p) = 1 − λ_{2, 3}(n, p), λ_{3, 3}(n, p) = 1 − λ_{3, 4}(n, p), λ_{4, 4}(n, p) = 1 − λ_{4, 1}(n, p) and λ_{1, 1}(n, p) = 1, leaving us to define λ_{2, 3}(n, p), λ_{3, 4}(n, p) and λ_{4, 1}(n, p). The resulting transition matrix is depicted in Table 3. It follows from the formulation of our model that an uninterrupted spell of proximity days between two public holidays will be deemed post-, rather than pre-, holiday days. This is to avoid difficulties in identifying a time of transition from the post- to the pre-holiday state.

Expert judgement from engineers at NGN suggested that a proximity effect, if it was felt, was likely to last for 1 or 2 days on either side of a public holiday. Longer periods of protracted holiday behaviour were thought to be very unlikely. If a pair of public holidays was separated by a small number of ordinary days, transitions into the proximity state, and transitions sustaining it, were expected to be more likely. To build the first of these ideas into our model, we allow the probability of transition into the pre-holiday state, λ_{4, 1}(n, p), to decay quickly to 0 as n becomes large by modelling the logit of the transition probability as

$$\text{logit}\{{\lambda }_{4,1}(n,p)\}={\nu }_{4,1,1}+{\nu }_{4,1,2}\frac{{(n-1)}^{1/2}}{10}.$$

(1)

We expect that ν_{4, 1, 2} < 0. Here the offset allows the parameter ν_{4, 1, 1} to be interpreted as the logit probability of transition into the pre-holiday state when the following day is a holiday. Taking the square root of n − 1 helps to prevent the prior variance of λ_{4, 1}(n, p) from becoming too large as n grows. In turn, this prevents an implausible U-shaped prior for any feasible n. Given the values that are taken by n, the division of (n − 1)^1/2 by 10 enables the coefficient ν_{4, 1, 2} to be interpreted on a scale that is comparable with that of a binary predictor, which will be introduced when modelling the other two transition probabilities, λ_{3, 4}(n, p) and λ_{2, 3}(n, p).

For the probability of transition out of the post-holiday state, one could assume that λ_{3, 4}(n, p) = λ_{3, 4} for all $(n,p)\in {\mathbb{N}}_0^2\backslash \{(0,0)\}$ and appeal to the property of (homogeneous) first-order Markov chains that the sojourn time in state 3 would have a geometric distribution with mean ${\lambda }_{3,4}^{-1}$⁠. However, a geometric distribution of the sojourn time in the post-holiday state, having only one parameter, is not sufficiently flexible to be reconciled with the expert judgement of NGN engineers. We therefore use the information that is encoded in the explanatory variable p, the number of days since the previous holiday, and model the logit of the transition probability as

$$\text{logit}\{{\lambda }_{3,4}(n,p)\}={\nu }_{3,4,1}+{\nu }_{3,4,2}\frac{{(p-2)}^{1/2}}{10}+{\nu }_{3,4,3}\mathbb{I}(n=1)$$

(2)

in which $\mathbb{I}(n=k)$ takes the value 1 if n = k and 0 otherwise. Clearly the transition probability λ_{3, 4}(n, p) will approach 1 as p becomes large if ν_{3, 4, 2} > 0. The offset in the square-root term enables the parameter ν_{3, 4, 1} to be interpreted as the logit probability of transition out of the post-holiday state after a stay of 1 day. In the UK, there is at least one weekend per year—the Easter weekend—with public holidays on the neighbouring Friday and Monday. We include the indicator term to allow the transition probability to differ over such weekends, in accordance with the views of industry experts at NGN. An analogous term is included in the logit probability of transition into the post-holiday state, which we model as

$$\text{logit}\{{\lambda }_{2,3}(n,p)\}={\nu }_{2,3,1}+{\nu }_{2,3,2}\mathbb{I}(n=2).$$

(3)

We have used a logit link function. As with any choice of a functional form in a regression model, there are other possibilities. However, in our application we found that the posterior means of the transition probabilities tended either to be very small or to lie in a fairly central range, roughly between 0.25 and 0.65, over which the probit link function, for example, can match the logit function closely and even an asymmetric link, such as the complementary log–log-function, can match quite closely. The greatest discrepancy between these three functions would then occur in the case of the complementary log–log-function deviating from the other two for larger probabilities, which tended not to occur. Nevertheless, there would be scope for future examination of this choice, given sufficient data.

3.2 Modelling the conditional demand for gas

Analysis involving an earlier version of the model, where observations were assumed to be conditionally independent given the states, revealed substantial auto-correlation between residuals, causing difficulties in identifying the three unknown states. Tunnicliffe-Wilson et al. (2016), section 3.10.4, also noted a carry-over effect in gas demand data. We therefore model the logarithm Y_t of the demand for gas, conditionally on the state on the current and previous days, as a first-order vector auto-regression, or VAR(1), with density $p({\mathbf{y}}_t|{\mathbf{y}}_{t-1},S_{t-1}=s_{t-1},S_t=s_t,n_t,p_t,r_t,{\mathbf{w}}_t)$ for t = 2,…, T and, at time t = 1,

$$p({\mathbf{y}}_1|S_0=j,S_1=s_1,n_1,p_1,r_1,{\mathbf{w}}_1)=p({\mathbf{y}}_1|S_1=s_1,n_1,p_1,r_1,{\mathbf{w}}_1).$$

(4)

We constrain our conditional model for gas demand to be conditionally stationary, i.e. we can regard the vectors Y_t as being the results of transforming values from a stationary VAR(1) process. The transformation applied on day t depends on the day of the year, day of the week, composite weather variable and state. Given the relatively short timescales of interest, such a stationarity assumption is reasonable and, in fact, non-stationarity of a non-negligible degree seems implausible in this application.

The conditional model for $({\mathbf{Y}}_t|{\mathbf{Y}}_{t-1}={\mathbf{y}}_{t-1},S_{t-1}=s_{t-1},S_t=s_t)$ is given by

$${\mathbf{Y}}_t-{\mathit{\mu }}_t=\mathrm{\Psi }({\mathbf{y}}_{t-1}-{\mathit{\mu }}_{t-1})+{\mathit{ϵ}}_t,{\mathit{ϵ}}_t\sim N_2(\mathbf{0},{\mathrm{\Omega }}_t^{-1}),$$

(5)

for t = 2,…, T, whereas, for the initial distribution of $({\mathbf{Y}}_1|S_1=s_1)$ in equation (4), we write

$${\mathbf{Y}}_1={\mathit{\mu }}_1+{\mathit{ϵ}}_1,{\mathit{ϵ}}_1\sim N_2\{\mathbf{0},V(\mathrm{\Psi },{\mathrm{\Omega }}_1)\}.$$

(6)

In both cases the dependence on the states comes through both the time-dependent mean μ_t and the precision matrix Ω_t. The terms ε₁,…,ε_T form a sequence of independent bivariate normal random vectors with zero mean, $\mathrm{\Psi }\in {\mathcal{S}}_{\psi }$ is a 2 × 2 real-valued matrix, with elements Ψ_{j, k}, constrained to satisfy the stationarity condition of a bivariate AR(1) process, and each Ω_t is a 2 × 2 symmetric, positive definite matrix. The stationarity region ${\mathcal{S}}_{\psi }$ is defined as the subset of 2 × 2 matrices with real-valued entries whose eigenvalues are less than 1 in modulus (Tunnicliffe-Wilson et al., 2016). The term V(Ψ,Ω₁) is the stationary variance of the mean-centred errors {Y_t−μ_t:t = 1,2,…} that would prevail if the precision matrix remained equal to Ω₁ at all future times, i.e. the matrix V which solves the equation $V=\mathrm{\Psi }V{\mathrm{\Psi }}^{\prime }+{\mathrm{\Omega }}_1^{-1}$⁠. A solution, which is symmetric and positive definite, is guaranteed by our assumptions on the support of Ψ and Ω₁. Specification of a well-behaved prior distribution over the stationarity region ${\mathcal{S}}_{\psi }$ is very challenging because of its complex geometric constraints. Fortunately, the problem simplifies greatly if we assume that Ψ_{1, 1} = Ψ_{2, 2} = Ψ_on and Ψ_{1, 2} = Ψ_{2, 1} = Ψ_off, in which case the necessary and sufficient conditions for stationarity are that $|{\mathrm{\Psi }}_{\text{on}}+{\mathrm{\Psi }}_{\text{off}}|<1$ and $|{\mathrm{\Psi }}_{\text{on}}-{\mathrm{\Psi }}_{\text{off}}|<1$⁠. We therefore adopt this simplification to the model.

The time-dependent mean μ_t in expressions (5) and (6) includes terms to allow for the influence of the CWV and various seasonal- and calendar-related effects. As we are working on the logarithmic scale, each has a multiplicative effect on the demand for gas. Conditionally on the state, S_t = s_t, the mean for LDZ j is given by

$${\mu }_{t,j}={\alpha }_j+B_{t,j,s_t}{\beta }_{j,r_t}+{\mathrm{\Gamma }}_{t,j}+{\mathrm{\Delta }}_{t,j}+({\zeta }_{j,1}+{\zeta }_{j,2}{w}_{t,j}){\tilde{w}}_{t,j}$$

(7)

for j = 1,2, where α_j provides an intercept. Because of interactions with the temperature and differences in consumer habits, the effect of a public holiday may differ according to whether it falls over the Christmas period, at Easter or over the summer. We therefore allow the three types of public holiday, which were defined earlier in Table 2, to have different effects on μ_{t, j}, represented by β_j,1, β_j,2 and β_j,3. The mean depends on the state on day t through the term $B_{t,j,s_t}$ which controls whether or not a holiday effect is included. It is defined by

$$\begin{array}{cc}& \hfill B_{t,j,1}={\rho }_{\beta ,j}^{n_t},\hfill \\ & \hfill B_{t,j,2}=1,\hfill \\ & \hfill B_{t,j,3}={\rho }_{\beta ,j}^{\text{min}(n_t,p_t)},\hfill \\ & \hfill B_{t,j,4}=0,\hfill \end{array}$$

where ${\rho }_{\beta ,j}\in (0,1)$ so that the holiday effect ${\beta }_{j,r_t}$ is scaled down on proximity days (states 1 and 3) by a factor which decays to 0 with increasing separation. Since the other terms in equation (7) do not depend on the state on day t, the mean μ_{t, j} in the pre- or post-holiday state is a weighted average of the mean in the normal and holiday states. The weights depend on the number of days separating day t from the next or closest public holiday.

In equation (7), the covariate ${\tilde{w}}_{t,j}$ is defined by ${\tilde{w}}_{t,j}=w_{t,j}-m_{d(t),j}$⁠, where m_{d(t), j} is a smoothed average for the CWV in region j on day t where d(t) ∈ {1,2,…,366}. We include both the mean-centred CWV, ${\tilde{w}}_{t,j}$⁠, and its interaction with the raw CWV, w_{t, j}, to allow the effect of above or below average temperatures to differ according to absolute weather conditions. The structure of this term is supported by the on-line supplementary Fig. S1 which shows the effect of w_{t, j} on the relationship between log(gas demand) and the mean-centred CWV.

The term Δ_{t, j} is constructed to give a day of the week effect, whereas Γ_{t, j} provides a seasonal component, after allowing for the CWV. Both terms are composed by using Fourier series with

$${\mathrm{\Delta }}_{t,j}=\sum _{k=1}^3\left\{{\delta }_{j,1,k}\cos \right(\frac{2\pi kt}{7})+{\delta }_{j,2,k}\sin (\frac{2\pi kt}{7}\left)\right\}$$

(8)

and

$${\mathrm{\Gamma }}_{t,j}=\sum _{k=1}^{K_{\gamma }}\left\{{\gamma }_{j,1,k}\cos \right(\frac{2\pi kt}{365.25})+{\gamma }_{j,2,k}\sin (\frac{2\pi kt}{365.25}\left)\right\}.$$

(9)

The six unconstrained parameters ${\mathit{\delta }}_j={({\delta }_{j,1,1},\dots ,{\delta }_{j,2,3})}^{\prime }\in {\mathbb{R}}^6$ in equation (8) provide fixed effects for each day of the week which, by construction, sum to 0. The advantages of this parameterization are that the six elements of δ_j are unconstrained and, by treating them as exchangeable a priori, we can induce a prior for the seven fixed effects which is symmetric with respect to the day of the week. This parameterization and prior do not appear to have been used before in the gas demand literature. The Fourier series for the seasonal term in equation (9) is truncated at a modest number K_γ of harmonics beyond which we assume that their contribution is negligible. We note that an alternative representation of Γ_{t, j} using cyclic cubic regression splines proved to be less successful at explaining the seasonal variation when using a small number of degrees of freedom. The choice of K_γ will be discussed further in Section 6.

For the precision matrix Ω_t in equations (5) and (6) we adopt an approach which we believe is novel in the gas demand literature. We allow Ω_t to depend on S_t and hence to vary over time, with different values for the holiday and normal states when S_t = 2 and S_t = 4. Adopting an approach that is similar to that taken for the mean μ_t, we model the precision matrix on proximity days (states 1 and 3) by interpolating between these two values, with weights that depend on the number of days to or from the neighbouring public holiday. This is most natural when working in ${\mathbb{R}}^3$⁠, rather than the constrained space of 2 × 2 covariance matrices. We therefore reparameterize the precision matrix Ω_t in terms of its square-root-free Cholesky decomposition ${\mathrm{\Omega }}_t=T_{{\mathrm{\Omega }}_t}^{\prime }{D}_{{\mathrm{\Omega }}_t}^{-1}{T}_{{\mathrm{\Omega }}_t}$ (Pourahmadi, 1999), in which $T_{{\mathrm{\Omega }}_t}$ is a unit lower triangular matrix with (2,1)-element − ϕ_t and $D_{{\mathrm{\Omega }}_t}=\text{diag}({\tau }_{t,1}^{-1},{\tau }_{t,2}^{-1})$⁠. The new parameters have a convenient interpretation in terms of the auto-regression of ε_t,2 on ε_t,1, with ${\varphi }_t\in \mathbb{R}$ representing the auto-regressive coefficient and τ_t,2 > 0 the associated conditional precision. The parameter τ_t,1 > 0 represents the marginal precision of ε_t,1. For the coefficient, ω_t,1 = ϕ_t, and the log-precisions, ω_t,2 = ln (τ_t,1) and ω_t,3 = ln (τ_t,2), we adopt linear models with, for i = 1,2,3,

$${\omega }_{t,i}={\eta }_i+{\mathrm{\Theta }}_{t,i,s_t}{\theta }_i+K_{t,i},$$

(10)

in which η_i is an intercept and θ_i allows for the effect of a public holiday. Again, we assume that days before and after a holiday have the same effect and define

$$\begin{array}{cc}& \hfill {\mathrm{\Theta }}_{t,i,1}={\rho }_{\theta }^{n_t},\hfill \\ & \hfill {\mathrm{\Theta }}_{t,i,2}=1,\hfill \\ & \hfill {\mathrm{\Theta }}_{t,i,3}={\rho }_{\theta }^{min(n_t,p_t)},\hfill \\ & \hfill {\mathrm{\Theta }}_{t,i,4}=0,\hfill \end{array}$$

where ${\rho }_{\theta }\in (0,1)$⁠, which scales down the holiday effect θ_i on proximity days. The term K_{t, i} in equation (10) provides a seasonal component. A seasonal component was deemed necessary after analysis with an earlier version of the model, in which K_{t, i} had been omitted, revealed seasonal variation in the residual variance, which caused confounding between the state allocation and low frequency seasonal change. As with the seasonal component Γ_{t, j} in the time-varying mean, we use a truncated Fourier series to represent K_{t, i}, taking

$$K_{t,i}=\sum _{k=1}^{K_{\kappa }}\left\{{\kappa }_{i,1,k}\cos \right(\frac{2\pi kt}{365.25})+{\kappa }_{i,2,k}\sin (\frac{2\pi kt}{365.25}\left)\right\}.$$

(11)

The choice of truncation point K_κ will be discussed further in Section 6.

4 Prior distribution

Denoting the unknown parameters of the Markov model for state evolution by Λ and the parameters of the conditional model for demand by Π, we adopt a prior distribution in which Λ and Π are independent. The two independent components of our prior are described in the sections which follow.

4.1 Transition probabilities

The parameters Λ=(ν_{4, 1, 1}, ν_{4, 1, 2}, ν_{3, 4, 1}, ν_{3, 4, 2}, ν_{3, 4, 3}, ν_{2, 3, 1}, ν_{2, 3, 2})^′ comprise the linear coefficients from the logit probabilities (1)–(3). In the absence of a justification for correlation in either direction in our prior beliefs about the values of these parameters, we adopt a prior with independence between the ν_{j, k, i} and take

$${\nu }_{j,k,i}\sim N(m_{j,k,i},v_{j,k,i}).$$

(12)

The symmetric logit normal distribution, created by assigning a normal N(0, v) distribution to the logit transformation of a probability-valued random variable, becomes bimodal when v > 1. U-shaped priors for the transition probabilities λ_{j, k}(n, p) were not consistent with our prior beliefs and so we moderated our choice of prior variances v_{j, k, i} to avoid bimodality. To avoid the complicated likelihood leading to implausible sets of parameter values having undue weight in the posterior, a process of trial and improvement, guided by the interpretation of the ν_{j, k, i}, was then used to choose the hyperparameters. For a recent discussion of such issues see, for example, Gelman et al. (2017).

4.2 Parameters of the conditional demand model

For region j = 1,2, let β_j = (β_j,1,…, β_j,3)^′, ζ_j = (ζ_j,1, ζ_j,2)^′, ${\mathit{\gamma }}_j={({\gamma }_{j,1,1},\dots ,{\gamma }_{j,1,K_{\gamma }},{\gamma }_{j,2,1},\dots ,{\gamma }_{j,2,K_{\gamma }})}^{\prime }$ and, similarly, δ_j = (δ_j,1,1,…, δ_j,2,3)^′. Define $\mathit{\beta }={({\mathit{\beta }}_1^{\prime },{\mathit{\beta }}_2^{\prime })}^{\prime }$⁠, with corresponding concatenated vectors for the other region-specific parameters. In the model for the time-varying precision matrix, let η=(η₁, η₂, η₃)^′, θ=(θ₁, θ₂, θ₃)^′, ${\mathit{\kappa }}_i={({\kappa }_{i,1,1},\dots ,{\kappa }_{i,2,K_{\kappa }})}^{\prime }$ and $\mathit{\kappa }={({\mathit{\kappa }}_1^{\prime },{\mathit{\kappa }}_2^{\prime },{\mathit{\kappa }}_3^{\prime })}^{\prime }$⁠. Finally, let ρ=(ρ_β,1, ρ_β,2, ρ_θ)^′ for the parameters governing the rate of decay of the holiday effects in the time-varying mean and precision matrix. The model parameters then comprise Π={Ψ,α,β,γ,δ,ρ,ζ,η,θ,κ}. We impose prior independence between these parameter blocks.

As discussed in Section 3.2 we assume that the auto-regressive coefficient matrix Ψ is composed of a common diagonal element Ψ_on and a common off-diagonal element Ψ_off so that the stationarity condition reduces to $|{\mathrm{\Psi }}_{\text{on}}+{\mathrm{\Psi }}_{\text{off}}|<1$ and $|{\mathrm{\Psi }}_{\text{on}}-{\mathrm{\Psi }}_{\text{off}}|<1$⁠. Our novel approach makes it easier to impose this constraint. For this, we first define χ₁ = Ψ_on + Ψ_off and χ₂ = Ψ_on−Ψ_off and then reparameterize Ψ in terms of new parameters ${\xi }_i=({\chi }_i+1)/2\in (0,1)$ for i = 1,2. To construct our prior, we make ξ₁ and ξ₂ independent and give them beta distributions, ξ_i∼beta(a_{ξ, i}, b_{ξ, i}).

For the parameters in the time-varying mean μ_t governing the intercept, influence of the CWV, and the seasonal, day of the week and public holiday effects, we adopt hierarchical priors which allow information to be shared between the NO (j = 1) and NE (j = 2) LDZs. For the intercept we take ${\alpha }_j|{\mu }_{\alpha }{\sim }^{\text{IID}}N\{{\mu }_{\alpha },(1-r_{\alpha })v_{\alpha }\}$⁠, for j = 1,2, with μ_α∼N(m_α, r_γv_α), in which the correlation between sites is r_α, the marginal variance is v_α and the marginal mean is m_α. Similarly, independently for the coefficients of the mean-centred CWV (i = 1) and its interaction with the raw value (i = 2), we choose ${\zeta }_{j,i}|{\mu }_{\alpha }{\sim }^{\text{IID}}N\{{\mu }_{\zeta ,i},(1-r_{\zeta ,i})v_{\zeta ,i}\}$⁠, for j = 1,2, with μ_{ζ, i}∼N(m_{ζ, i}, r_{ζ, i}v_{ζ, i}). For the day of the week effects, independently for m = 1,2 and k = 1,2,3 we choose ${\delta }_{j,m,k}|{\mu }_{\delta ,m,k}{\sim }^{\text{IID}}N\{{\mu }_{\delta ,m,k},(1-r_{\delta })v_{\delta }\}$⁠, for j = 1,2, with μ_{δ, m, k}∼N(0, r_δv_δ). For the seasonal effects, independently for m = 1,2 and k = 1,…, K_γ, we take ${\gamma }_{j,m,k}|{\mu }_{\gamma ,m,k}{\sim }^{\text{IID}}N\{{\mu }_{\gamma ,m,k},(1-r_{\gamma })v_{\gamma ,k}\}$⁠, for j = 1,2, with μ_{γ, m, k}∼N(0, r_γv_{γ, k}), in which $v_{\gamma ,1}⩾v_{\gamma ,2}⩾\dots ⩾v_{\gamma ,K_{\gamma }}$⁠. We note that the assignment of independent, zero-mean normal distributions N(0, v_{γ, k}) to a pair of Fourier coefficients γ_{j,1, k} and γ_{j,2, k} is equivalent to assigning a uniform distribution to the phase of the kth harmonic and, independently, a Rayleigh distribution to its amplitude, with scale parameter √v_{γ, k}. Our prior therefore conveys the idea that the size of the seasonal harmonics will decay as their frequency increases.

For the public holiday effects, we adopt a prior of the form ${\mathit{\beta }}_j|{\mathit{\mu }}_{\beta }{\sim }^{\text{IID}}{N}_3\{{\mathit{\mu }}_{\beta },(1-r_{\beta }){\mathbf{V}}_{\beta }\}$⁠, For j = 1,2, with μ_β∼N₃(0, r_βV_β). Here V_β is a 3 × 3 compound symmetric matrix with non-zero off-diagonal elements, allowing positive correlation between the effects of each type of public holiday so that information can be pooled across types. This offers a compromise between an impractical assumption that the β_{j, i} are independent for each site j, which would limit the information that is available for inference, and an inflexible assumption that they are all equal.

The factor ρ_{β, j} acts on proximity days to scale down the holiday effect in the time-varying mean μ_{t, j} for LDZ j, j = 1,2. Similarly, the factor ρ_θ scales down the holiday effect in our reparameterization of the time-varying precision matrix, Ω_t, of the errors ε_t. The factors ρ_β,1 and ρ_β,2 in the means for the NO and NE LDZs have the same function and we expect them to be very similar. Although we would also expect the ρ_{β, j} to be informative about ρ_θ, we believe that they will carry less information because they relate to the rate of decay of parameters which play different roles. We therefore construct an asymmetric hierarchical prior by taking ${\tilde{\rho }}_{·}=\text{logit}({\rho }_{·})=\ln \{{\rho }_{·}/(1-{\rho }_{·})\}$⁠, and then choosing ${\tilde{\rho }}_{\beta ,j}|{\tilde{\rho }}_{\beta }{\sim }^{\text{IID}}N\{{\tilde{\rho }}_{\beta },(1-r_{\tilde{\rho },1})v_{\tilde{\rho }}\}$ in the means for LDZs j = 1,2 and then ${\tilde{\rho }}_{\beta }|{\mu }_{\tilde{\rho }}\sim N\{{\mu }_{\tilde{\rho }},(1-r_{\tilde{\rho },2})r_{\tilde{\rho },1}{v}_{\tilde{\rho }}\}$ and ${\tilde{\rho }}_{\theta }|{\mu }_{\tilde{\rho }}\sim N\{{\mu }_{\tilde{\rho }},(1-r_{\tilde{\rho },2})r_{\tilde{\rho },1}{v}_{\tilde{\rho }}\}$ and finally ${\mu }_{\tilde{\rho }}\sim N(m_{\tilde{\rho }},r_{\tilde{\rho },1}{r}_{\tilde{\rho },2}{v}_{\tilde{\rho }})$⁠. The marginal prior moments for the logit parameters are $E({\tilde{\rho }}_{\beta ,j})=E({\tilde{\rho }}_{\theta })=m_{\tilde{\rho }}$⁠, $\text{var}({\tilde{\rho }}_{\beta ,j})=v_{\tilde{\rho }}$⁠, $\text{corr}({\tilde{\rho }}_{\beta ,j},{\tilde{\rho }}_{\theta })=\surd r_{\tilde{\rho },1}{r}_{\tilde{\rho },2}$⁠, for j = 1,2, and $\text{var}({\tilde{\rho }}_{\theta })=r_{\tilde{\rho },1}{v}_{\tilde{\rho }}$ and $\text{corr}({\tilde{\rho }}_{\beta ,1},{\tilde{\rho }}_{\beta ,2})=r_{\tilde{\rho },1}$ where $r_{\tilde{\rho },1},r_{\tilde{\rho },2}\in (0,1)$⁠. We can therefore choose the correlation between the ${\tilde{\rho }}_{\beta ,j}$ to be greater than their correlation with ${\tilde{\rho }}_{\theta }$⁠.

Finally, the intercept and holiday effects in our reparameterized time-varying precision matrix are given independent normal distributions, taking η_i∼N(m_{η, i}, v_{η, i}) and θ_i∼N(0, v_{θ, i}) for i = 1,2,3. Similarly, the seasonal effects are given independent normal priors, with κ_{i, m, k}∼N(0, v_{κ, i, k}), m = 1,2, k = 1,…, K_κ, for i = 1,2,3.

Our choices for the hyperparameters in each component of our prior distribution are detailed in Section 6.1.

5 Posterior inference

The posterior distribution of the model parameters follows from Bayes theorem as

$$\pi (\mathrm{\Pi },\mathrm{\Lambda }|\mathbf{y},\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})\propto p(\mathbf{y}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})\pi (\mathrm{\Pi })\pi (\mathrm{\Lambda })$$

in which y represents the complete time series, here y=y_1:T, where the general notation x_i:j for i<j indicates a sequence (x_i, x_i+1,…, x_{j − 1}, x_j). The term p(y|Π,Λ,n,p,r,w) is the observed data likelihood given by

$$p(\mathbf{y}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})=\sum _{\mathbf{s}}p(\mathbf{y}|\mathbf{s},\mathrm{\Pi },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})p(\mathbf{s}|\mathrm{\Lambda },\mathbf{n},\mathbf{p})$$

in which the sum is taken over all possible sequences of states, s=s_0:T. It can be calculated efficiently by using a filtering algorithm which computes $\text{Pr}(S_t=k,{\mathbf{y}}_{1:t}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})$⁠, $k\in {\mathcal{S}}_s$⁠, recursively for t = 0,…, T and finally

$$p(\mathbf{y}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})=\sum _{k=1}^4\text{Pr}(S_T=k,{\mathbf{y}}_{1:T}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w}).$$

The lagged dependence of Y_t+1 on S_t in the directed acyclic graph in Fig. 1 complicates recursive algorithms for quantifying posterior uncertainty in the hidden states. However, the model can be reformulated to simplify its conditional independence structure by defining an augmented state on day t by ${\tilde{S}}_t={(S_{t-1},S_t)}^{\prime }$ for t = 1,…, T. The demand Y_t is then conditionally independent of the augmented state on the previous day ${\tilde{S}}_{t-1}$ given the augmented state on the current day ${\tilde{S}}_t$ and the previous observation Y_{t − 1}. As some transitions in our original state process have zero probability, there are only 11, rather than 4²=16, possible values for the augmented state ${\tilde{S}}_t$ and we denote this set of values by ${\mathcal{S}}_{\tilde{s}}$⁠. The mapping from ${\mathcal{S}}_s\times {\mathcal{S}}_s$ to ${\mathcal{S}}_{\tilde{s}}$ is detailed in Table S1 of the on-line supplementary materials. The transition matrix for the new hidden process is sparse and can be represented in terms of the original transition probabilities λ_{j, k}(n, p), as indicated in supplementary Table S2. The initial distribution $\text{Pr}({\tilde{S}}_1|\mathrm{\Lambda },{\mathbf{p}}_{0:1},{\mathbf{n}}_{0:1})$ can similarly be formed by taking an appropriate product of terms, l_j(n₀, p₀) and λ_{j, k}(n₁, p₁). Filtering algorithms for hidden Markov models of this simple form are standard. See, for example, chapter 2 of MacDonald and Zucchini (1997). Our algorithm is given in section S2.2 of the supplementary materials.

The posterior distribution π(Π,Λ|y,n,p,r,w) cannot be evaluated in closed form and so we build a numerical approximation using a Hamiltonian Monte Carlo (HMC) sampling scheme. We note that the likelihood p(y|Π,Λ,n,p,r,w) is not symmetric with respect to the labels of the hidden states. As a result the state-specific parameters are all identifiable in the posterior, and so methods to force an identified sample, such as those discussed in Stephens (2000), are not required.

5.1 Approximating the posterior for the model parameters

HMC sampling is a special case of the Metropolis algorithm; for example, see Neal (2011) or Betancourt (2017) for an introduction. It relies on the introduction of auxiliary variables that are interpreted as the momentum of a particle whose position in space is represented by the parameter values. This enables efficient proposals to be generated by exploiting Hamiltonian dynamics and modelling the movement of the sampler around the joint posterior of the momentum and position variables as the motion of the particle through an unbounded, frictionless space. We use rstan (Stan Development Team, 2016), the R interface to the Stan software (Carpenter et al., 2017), to implement the HMC algorithm. Stan requires users to write a program in the probabilistic Stan modelling language, the role of which is to provide instructions for computing the logarithm of the kernel of the posterior density function. The Stan software then automatically tunes and runs a Markov chain simulation to sample from the resulting posterior.

5.2 Approximating the posterior for the hidden states

Although the hidden states S_t are not sampled as part of the HMC scheme, their marginal posteriors can be approximated by Rao–Blackwellization through

$$\widehat{\text{Pr}}(S_t=k|\mathbf{y},\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})=\frac{1}{M}\sum _{i=1}^M\text{Pr}(S_t=k|\mathbf{y},{\mathrm{\Pi }}^{[i]},{\mathrm{\Lambda }}^{[i]},\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w}),k\in {\mathcal{S}}_s,$$

(13)

for t = 0,…, T where Π^[i] and Λ^[i] denote the ith posterior samples of the model parameters Π and Λ for i = 1,…, M. The full sample smoothed probabilities $\text{Pr}({\tilde{S}}_t=k|\mathbf{y},\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})$ can be computed in a backward recursion, starting at time t = T, using the forward probabilities $\text{Pr}({\tilde{S}}_t=k,{\mathbf{y}}_{1:t}|\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})$ calculated when computing the observed data likelihood. We can then marginalize over S_{t − 1} to compute the probabilities in our original four-state model $\text{Pr}(S_t=k|\mathbf{y},\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})$⁠, $k\in {\mathcal{S}}_s$⁠, for use in expression (13). The complete algorithm is given in section S2.3 of the on-line supplementary materials.

6 Application to daily demand data

The daily gas consumption data from NGN were introduced in Section 2. Using the algorithm that was described in Section 5, we fitted our hierarchical model to data from each of the three NDM load bands. When including a large number of Fourier components in the model for mean log-demand, the amplitude ${({\gamma }_{j,1,k}^2+{\gamma }_{j,2,k}^2)}^{1/2}$ of the kth harmonic in expression (9) for each site j was negligible when k > 6 and so we truncated the series at K_γ = 6. Applying similar reasoning, we truncated the Fourier series in expression (11) for the seasonally varying precision at K_κ = 12.

6.1 Prior specification

The structure of the prior distribution for the model parameters was outlined in Section 4. It began with an assumption of independence between the parameters Λ of the marginal model for the hidden states and the parameters Π of the conditional model for gas demand given the states. The moments in the prior for Π were chosen according to the specification in section S3.1 of the on-line supplementary materials and the densities were generally very flat.

Completing our prior for Λ requires choices of means and variances in the normal distributions (12) for the logit coefficients. Our general procedure was described in Section 4.1. By several iterations of calculation and inspection we arrived at the specification in section S3.2 of the supplementary materials which yielded the distribution for the states shown in Fig. 5 later in Section 6.3.2 over a representative year. The pointwise prior mode for the entire state sequence is displayed in supplementary Fig. S2.

6.2 Hamiltonian Monte Carlo implementation

The Stan program representing our hierarchical model is available from

https://academic.oup.com/jrsssa/issue.

For the data from each load band, using the rstan interface to the Stan software, we ran four HMC chains initialized at different starting points for 10000 iterations, half of which were discarded as burn-in. It was convenient to compute the posterior samples for the smoothed probabilities $\text{Pr}(S_t=l|\mathbf{y},\mathrm{\Pi },\mathrm{\Lambda },\mathbf{n},\mathbf{p},\mathbf{r},\mathbf{w})$ in expression (13) on line. To reduce the associated storage overheads, we therefore thinned the output to retain every 10th posterior draw. The usual graphical and numerical diagnostics gave no evidence of any lack of convergence and the effective sample size was at least 1565 for every parameter.

6.3 Posterior inference

6.3.1 Parameter inference

Fig. 2 shows the fixed effects β_{j, k} in the mean μ_{t, j} (7) for each type of public holiday in each LDZ. As expected, there is evidence that the mean takes larger values on public holidays in load band 1, which represents domestic customers, and smaller values in load bands 2 and 3, which broadly represent commercial and small industrial customers, whose business activity typically stops on holidays. The effect seems to be the largest among industrial customers (load band 3) and smallest among domestic users. In general, the effect of a public holiday is less pronounced during the public holidays around Easter. For each load band, if we compare the effects across the two LDZs, we can see that the holiday effects in load band 1 are nearly identical, which suggests that any differences in load bands 2 and 3 arise due to differences in the mix of customers that they represent. Different types of, for example, commercial, educational and industrial customers may have different patterns of closure. In fact some, such as shops, may remain open on certain public holidays. The mix of such, especially industrial, customers will differ between the two LDZs. Moreover, there may be differences in customary closure patterns between, say, the West Yorkshire conurbation (NE) and Tyneside (NO). It is not surprising, therefore, that we observe differences between LDZs in the holidays effects in load bands 2 and 3, most notably the ‘Easter’ effect in load band 2 and the ‘other’ effect in load band 3. This is exemplified in the on-line supplementary Fig. S3 which plots the marginal posterior densities of the differences, β_{1, k} − β_{2, k}, in the effects of each type, k = 1,2,3, of public holiday for each load band.

Fig. 2

Marginal posterior densities for the fixed effects β_{j, k} in the mean μ_{t, j} for each type, k = 1,2,3, of public holiday in each LDZ, j = 1,2 for (a) Easter, (b) other and (c) Christmas, and (d) marginal posterior density for the parameter ρ_{β, j} governing the rate of decay of the holiday effect on proximity days in each LDZ, j = 1,2: , load band 1, NO; , load band 1, NE; , load band 2, NO; , load band 2, NE; , load band 3, NO; , load band 3, NE

Fig. 2(d) shows the marginal posterior density for the parameter ρ_{β, j} which governs the rate of decay of the holiday effect in the pre- and post-holiday states. The posteriors for each LDZ in load band 3 seem to support larger values than those for load bands 1 and 2, suggesting a slower rate of decay.

The corresponding plots for the fixed effects θ_i in the parameters ω_{t, i} (10) of the square-root-free Cholesky decomposition of the precision matrix Ω_t are shown in Fig. 3. For all three load bands, there is strong evidence of a positive effect on the auto-regressive coefficient ω_t,1 = ϕ_t, suggesting a stronger positive correlation between the residuals in the two LDZs on public holidays. For the logarithms of the marginal error precision ω_t,2 = ln (τ_t,1) and the conditional error precision ω_t,3 = ln (τ_t,2), there is little evidence of a public holiday effect, except for a negative relationship with the log-conditional-precision in load band 3. This suggests that public holidays are associated with higher residual variance for log(gas demand) by industrial customers.

Fig. 3

Marginal posterior densities for the fixed effects θ_i in the parameters ω_{t, i} of the square-root-free decomposition of the precision matrix Ω_t, ω_t,1 = ϕ_t, (a) ω_t,3 = ln (τ_t,2), the auto-regressive component, (b) ω_t,2 = ln (τ_t,1) the log-marginal precision, (c) and (d) the parameter ρ_θ governing the rate of decay of the holiday effect on proximity days: , load band 1; , load band 2; , load band 3

6.3.2 State inference

For load bands 1–3 respectively, Figs 4(a)–4(c) show the pointwise posterior mode for the state sequence S₁,…, S_T. Different patterns are evident across the different load bands. In load band 1, there is little evidence of a proximity effect, with the posterior probability of assignment to the pre- and post-holiday states, 1 and 3, being less than 0.5 in the neighbourhood of most public holidays. For load bands 2 and 3, Figs 4(b) and 4(c) provide clear evidence of a proximity effect around the Christmas period and over the Easter weekend, with additional evidence of a post-holiday effect after the late spring and summer public holidays in load band 2. The sojourn times in the proximity states are longest around Christmas, most notably before Christmas Day in load band 3 and after Boxing Day in load band 2.

Fig. 4

Pointwise posterior mode for the state sequence S_t for load bands (a) 1, (b) 2 and (c) 3: , state 1; , state 2; , state 3; , state 4

Clearly the pointwise posterior mode for the state sequence cannot give any indication of the uncertainty in the state allocation. Therefore, using load band 3 as an example, Fig. 5 shows the posterior distribution for the states S_t in a representative year (2015), with the prior distribution overlaid. There is a marked difference between the prior and posterior in the period around a public holiday. Although the prior treats all public holidays as if they are (nearly) exchangeable, there is strong evidence of differences between public holidays in the posterior. This shows the advantage of allowing the extent of proximity days to be unknown rather than periods of common, fixed duration around each public holiday. The corresponding plots for load bands 1 and 2 are given in Figs S4 and S5 of the on-line supplementary materials. Whereas load band 2 shows similar patterns to load band 3, there is very little difference between the prior and posterior for load band 1, though the prior does assign slightly more mass to the pre- and post-holiday states in the neighbourhood of public holidays. This suggests that there is little evidence about the presence, or otherwise, of a proximity effect in the data for load band 1.

Fig. 5

Prior and posterior for the state sequence S_t for load band 3 in an illustrative year (2015): , state 1, prior; , state 1, posterior; , state 2, prior; , state 2, posterior; , state 3, prior; , state 3, posterior; , state 4, prior; , state 4, posterior

6.3.3 Posterior predictive inference

Consider a simplified version of the model that was described in Section 3 which ignores the proximity effect and omits the pre- and post-holiday states (states 1 and 3) so that each day is classified either as a public holiday (state 2) or not (state 4). As discussed in Section 1, many models from the literature are structured in this way, allowing for public holidays but not any protracted effect on neighbouring days. To illustrate the benefit of incorporating the proximity effect, we use the framework of posterior predictive checks (Gelman et al., 2013) to compare inferences under our four-state non-homogeneous Markov model with those computed under the simplified two-state model. In this framework, the basic idea is to measure the extent to which a model captures some data summary of interest by comparing its posterior predictive distribution with the value that was observed. In our case the (approximate) posterior predictive distribution is computed numerically on the basis of a Markov chain Monte Carlo sample from the posterior of the model parameters by simulating replicated data sets in one-to-one correspondence with the posterior draws. If the model can capture adequately the aspect of the data summary of interest, the observed value will look plausible under its posterior predictive distribution.

On public holidays and days which are not in the neighbourhood of a public holiday, both the two-state and the four-state models seem to provide a good fit to the data for all load bands. For instance, in load band 3, for all the days during the observation period which were public holidays, Fig. 6 compares the posterior predictive distribution for log(gas demand) under each model with the values that were observed. Similarly, Fig. S6 in the on-line supplementary materials provides an analogous comparison for days which were 10 days away from a public holiday. In these plots, points lying above the unit diagonal indicate dates for which predicted log-demand exceeds observed log-demand, and conversely for points below the diagonal. Good model fit is indicated by the majority of points lying close to the unit diagonal, suggesting that most observations fall within the main body of their posterior predictive distribution. This is broadly true for both models in each LDZ, with only 1.15–6.67% of observations lying outside the central 95% of their associated posterior predictive distribution. Interestingly, in Fig. 6, most of the observations that are flagged as lying in the (upper) tail of their posterior predictive distribution arise from just two dates: the Royal wedding on April 29th, 2011, and the Diamond Jubilee of Queen Elizabeth II on May 6th, 2012. These were one-off public holidays, rather than annual events, and it appears that they were associated with less of a drop in the demand for gas in load band 3 than in other public holidays occurring around that time of the year. This may have been because the holidays were not observed by all industrial customers. In contrast, in load bands 1 and 2, the demand for natural gas on these public holidays was modelled well.

Fig. 6

For load band 3 and (a) LDZ NO and (b) LDZ NE, posterior predictive means versus observed log(gas demand) for each public holiday in the observation period (upper panels refer to the simple two-state model; lower panels to the four-state non-homogeneous Markov model): , , 2.5% and 97.5% points in the posterior predictive distributions; , observation inside the central 95% of the posterior predictive distribution; , observation outside the central 95% of the posterior predictive distribution

The advantage of allowing a proximity effect becomes clear when we consider days that are close to a public holiday. For example, Fig. 7 shows the corresponding posterior predictive checks for all the days which were 1 day from a public holiday in load band 3. For each LDZ, it appears that the two-state model systematically overestimates the demand for gas, with 9.63% and 10.67% of observations lying outside (and mostly below) the central 95% of their posterior predictive distributions in the NO and NE LDZs respectively. In contrast, the posterior predictive densities under the four-state model generally support smaller values and are more diffuse, reflecting the additional uncertainty in the demand for gas that typically accompanies the period around public holidays. Correspondingly, the observations appear more plausible under the posterior predictive distributions, with only 5.19% and 4.44% lying outside the central 95% of their posterior predictive distributions in the two LDZs. A similar effect, albeit less marked, is observed for days which are 2 or 3 days from a public holiday. Very similar conclusions can be drawn based on the analysis of the data from load band 2, plots for which are provided in Figs S7–S9 of the on-line supplementary materials. Not surprisingly, there is no notable difference between the two- and four-state models for load band 1, where our analysis suggested little evidence of a proximity effect. Supplementary Table S3 contains a summary of the results that are described in this section for all load bands.

Fig. 7

For load band 3 and (a) LDZ NO and (b) LDZ NE, posterior predictive means versus observed log(gas demand) for each day in the observation period which was 1 day from a public holiday (upper panels are based on the simple two-state model; lower panels on the four-state non-homogeneous Markov model): , , 2.5% and 97.5% points in the posterior predictive distributions; , observation inside the central 95% of the posterior predictive distribution; , observation outside the central 95% of the posterior predictive distribution

7 Discussion

The energy sector in the UK is changing to harness ecologically more sound technologies in the face of growing economic, environmental, societal and public health concerns. As a comparatively low cost source of energy, natural gas has an important role to play in this changing energy mix. Yet, clearly, its contribution is contingent on efficient operation of the gas industry. One of the ways in which operational efficiency and decision making, more generally, can be enhanced is through improved forecasting of the demand for gas.

We have presented a novel model for daily gas demand which allows for a protracted effect of public holidays which extends into neighbouring days. A key feature of the model is that we allow the existence, duration and location of the proximity days to be unknown. This is achieved by modelling the data by using a four-state non-homogeneous Markov model with cyclic dynamics, and states that represent pre- and post-holiday days, as well as (observable) public holidays and normal days. The dates of public holidays are always known in advance and we allow for this by making our model non-homogeneous, with transition probabilities that depend on the number of days to the next, and since the previous, public holiday. We establish an interpretable parameterization for the logit transition probabilities and illustrate how this can be used to assign a prior distribution for the states which represents a flexible compromise between a model which does not allow for a proximity effect and a model which fixes the dates over which a proximity effect is felt.

Conditionally on the states, we model the natural logarithm of gas demand over two large geographical regions by using a first-order vector auto-regression, which is conditionally stationary. A novel feature of this model is the assumption of a symmetric, though non-diagonal, auto-regressive coefficient matrix which greatly simplifies the geometry of the stationarity region, allowing it to be expressed as a unit square. Other features of our model which are, to the best of our knowledge, novel in the gas demand literature include the symmetric treatment of the day of the week effect and the incorporation of seasonal and holiday effects in the variance of the residuals.

We applied the model to data from NGN, which is responsible for the distribution of gas to most of the north-east of England, northern Cumbria and Yorkshire. An assessment of model fit by using posterior predictive checks highlighted the improvement that is gained by allowing for a proximity effect among commercial and small industrial customers. Similarly, the posterior distribution for the states revealed different patterns in the identification of proximity days between public holidays and across different groups of customers. This highlights the advantage of allowing uncertainty in the identification of proximity days compared with an inflexible approach that treats them as periods of common fixed duration around each public holiday. The results of a preliminary version of the model are already being used by NGN in their annual medium-term forecasting exercise, with plans to extend the approach to consider the demand for gas by larger industrial customers.

Possible enhancements to the model include allowing the transition probabilities to depend, not just on the proximity of the preceding and following holidays, but also on the type of holiday and seasonal or temperature effects. With the limited number of observed holidays in our data, inference for the additional parameters would be difficult but this remains a topic for future research.

Acknowledgements

The authors thank NGN for providing the data and for many helpful discussions and comments during the research project. We also extend our gratitude to Matt Linsley, Manager of the Industrial Statistics Research Unit at Newcastle University, for his important role in organizing and facilitating meetings with NGN. We are also grateful to two reviewers for comments which have led to improvements in the paper.