Analysis of residential gas consumption data in Luxembourg

List of abbreviations

Abbreviation	Description
ACF	Autocorrelation function
GWh	Gigawatt hour
ILR	Institut Luxembourgeois de Régulation
IQR	Interquartile Range
MAD	Mean absolute deviation
MSD	Mean squared deviation
Q1	First quartile
Q3	Third quartile
Var	Variance
sd	Standard deviation

1 Introduction

This document presents the analysis of residential gas consumption data in Luxembourg.

2 Data sources

The analysis hereafter is based on two data sources with the following characteristics.

Data source	Producer	URL	Format
Monthly residential gas consumption in Luxembourg	Institut Luxembourgeois de Régulation (ILR)	https://data.public.lu/en/datasets/donnees-statistiques-du-secteur-de-gaz-naturel-ilr/	.xlsx
Monthly mean air temperature in Luxembourg	MeteoLux	https://data.public.lu/en/datasets/inspire-annex-iii-meteorological-geographical-features-pointtimeseriesobservation-monthly-weather-measurements-at-luxembourg-findel-airport-1/	.gml

For conciseness, the details of data download, import and tidying are not discussed here. If you are interested, you can use the R code from GitHub to download, import, and tidy the data.

3 Data description

For our analysis, we are only interested in the monthly residential gas consumption in GWh. Note that the raw data has two other types of gas consumption (industrial consumption and production/co-generation).

We use tsibble objects as they are specifically designed for time series analysis. A tsibble is a special kind of tibble.

After downloading, importing, and tidying the data, the resulting tsibble, named df is shown below.

We observe that:

• There are 93 rows and 3 columns
• The interval (time step) is 1 month (month column)
• There are two time series:
  • Monthly residential gas consumption in GWh (consumption_gwh column)
  • Monthly average temperature in degrees Celcius (temp_celcius column)

Although the first 10 rows are displayed, the month column contains 8 years.

unique(year(df$month))

## [1] 2017 2018 2019 2020 2021 2022 2023 2024

Furthermore, we observe that the time series:

• starts at Jan-2017
• ends at Sep-2024
• has 0 gaps

Next, we perform an exploratory data analysis to understand the patterns in the data.

4 Exploratory data analysis

4.1 Desriptive statistics

The following descriptive statistics are computed for the monthly temperature in Celcius.

Statistic	Value
n	93.0
min	-1.6
q1	4.9
mean	10.8
median	10.0
q3	16.7
max	22.1

For the residential gas consumption in GWh, we compute a more comprehensive list of statistics.

consumption_gwh_stats <- df |>
  as_tibble() |>
  summarise(
    n = n(),
    sum = sum(consumption_gwh),
    min = min(consumption_gwh),
    q1 = quantile(consumption_gwh, 0.25),
    mean = mean(consumption_gwh),
    median = median(consumption_gwh),
    q3 = quantile(consumption_gwh, 0.75),
    max = max(consumption_gwh),
    var = var(consumption_gwh),
    sd = sd(consumption_gwh),
    skewness = skewness(consumption_gwh),
    kurtosis = kurtosis(consumption_gwh),
    IQR = IQR(consumption_gwh),
    MAD = mad(consumption_gwh, center = mean(consumption_gwh)),
    MSD = mean((consumption_gwh - mean(consumption_gwh))^2)
  ) |>
  pivot_longer(everything(),
    names_to = "Statistic",
    values_to = "Value"
  ) |>
  mutate(Value = round(Value, 2))


kbl(consumption_gwh_stats) |>
  kable_styling(bootstrap_options = "striped",
                full_width = FALSE,
                position = "center")

Statistic	Value
n	93.00
sum	21212.76
min	32.47
q1	172.14
mean	228.09
median	225.09
q3	309.17
max	416.42
var	10287.68
sd	101.43
skewness	-0.16
kurtosis	2.17
IQR	137.03
MAD	102.52
MSD	10177.06

We observe that:

• The total residential gas consumption was 21,212.76 GWh.
• The average monthly residential gas consumption is 228.09 GWh.

4.2 Histograms

Let us examine the distribution of the data.

4.3 Density plots

We observe that both distributions are bimodal. The temperature distribution has a more pronounced bimodal feature.

4.4 Boxplots

One simple yet effective method to detect outliers graphically is the boxplot. These are the boxplots for our two variables.

The lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The line between the hinges corresponds to the median. The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge. IQR is the distance between the first and third quartiles. The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data beyond the end of the whiskers are called “outliers” and are plotted individually.

We observe that both distributions do not contain any outliers.

4.5 Time plots

Perhaps the most useful data visualization when it comes to time series is the time plot. Here is the time plot for our data.

The time plot of the gas consumption (upper plot) immediately reveals some interesting features.

• Peaks happen in the winter.
• Low points happen in the summer.
• There is a decreasing trend.
• The seasonal pattern increases in size as the level of the series increases. This causes the variance to increase over time and it is also known as heteroskedasticity.

The spikes during the winter months are probably caused due to household gas heating. Similarly, the low points in the summer could be due to the reduced need for household gas heating.

The time plot of the mean monthly temperature (lower plot) exhibits the opposite pattern.

• Peaks happen in the summer
• Low points happen in the winter.
• There is a flat trend.
• The data appears to be homoskedastic.

Both time series have a strong seasonal pattern.

4.6 Seasonal plots

A seasonal plot is similar to a time plot except that the data are plotted against the individual month in which the data were observed.

In addition to the observations made by the time plot, we observe that:

• February and March 2018 have a higher gas consumption than the same months in other years.
• February 2024 has a lower gas consumption than the same month in other years.
• There are two “clusters” of curves. The first cluster has 2017 to 2021 curves and the second cluster has 2022 to 2024 curves. The curves of the second cluster are below than the curves of the first cluster, from April to November. This could be explained by the higher-than-usual residential gas prices observed from 2022 until

The higher-than-usual gas consumption in February and March 2018 is most likely explained by the low mean temperature during these months, as observed in the graph below.

4.7 Seasonal subseries plot

An alternative plot that emphasizes the seasonal patterns is where the data for each season are collected together in separate mini time plots. The blue horizontal lines indicate the means for each month. The x-axis indicates the year.

We observe that the pattern is similar between months, except from January.

We do not detect any clear pattern in this seasonal subseries graph.

Up until now, we have observed each time series individually. This is also called a univariate analysis. Next, we will look at how these time series behave together. This is also called a multivariate analysis.

4.8 Scatterplot

Let us visualize the data points of gas consumption and mean temperature together. The most common plot for this task is the scatterplot, as shown below. A linear fit is also plotted.

As expected, we observe a negative correlation between mean temperature and gas consumption. We also observe that the upper left corner (low temperature, high consumption) has data points closer to the linear regression line. In contrast, the lower right corner (high temperature, low consumption) has data points scattered away from the linear regression line.

The correlation coefficient measures the strength of the linear relationship between mean temperature and gas consumption. For this data, the correlation coefficient is -0.85.

As shown below, the correlation test indicates that the correlation coefficient is significantly different than zero.

cor.test(df$consumption_gwh, df$temp_celcius, method = "pearson",
         alternative = "two.sided", conf.level = 0.95)

## 
##  Pearson's product-moment correlation
## 
## data:  df$consumption_gwh and df$temp_celcius
## t = -15.48, df = 91, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8991106 -0.7834962
## sample estimates:
##        cor 
## -0.8513285

4.9 Lag plots

The figure below displays scatterplots where the horizontal axis shows lagged values (k) of the time series. Each graph shows
monthly gas consumption data points plotted against themselves at different lags.

df |>
  gg_lag(consumption_gwh, geom = "point", lags = 1:12) +
  labs(
    x = "lag(gas consumption in GWh, k)",
    y = "Gas consumption in GWh",
    color = "Month"
  ) +
  theme_minimal()

The colors indicate the month of the variable on the vertical axis. The relationship is strongly positive at lags 1 and 12, reflecting the strong seasonality in the data. The negative relationship seen for lag 6 occurs because peaks of gas consumption (in the winther) are plotted against troughs (in the summer).

4.10 Autocorrelation

Just as correlation measures the extent of a linear relationship between two variables, autocorrelation measures the linear relationship between lagged values of a time series.

There are several autocorrelation coefficients, corresponding to each panel in the lag plot above.

The autocorrelation coefficients make up the autocorrelation function or ACF.

df |>
  ACF(consumption_gwh, lag_max = 48) |>
  autoplot() +
  labs(
    title = "ACF plot: Monthly gas consumption in GWh",
    x = "Lag in months",
    y = "ACF"
  ) +
  theme_minimal()

In this graph:

• The autocorrelation coefficients at lag 1, 12, 24, etc. (denoted \(r_1\), \(r_{12}\), and \(r_{24}\) respectively) are high. This is due to the seasonal pattern in the data: the peaks tend to be 12 months apart and the low points tend to also be 12 months apart.
• The dashed blue lines indicate whether the correlations are significantly different from zero. If an autocorrelation coefficient lies outside the blue bound, it is significantly different from zero.

4.11 Box-Cox transformation

For more details, please refer to (Hyndman & Athanasopoulos, 2021).

lambda <- df |>
  features(consumption_gwh, features = guerrero) |>
  pull(lambda_guerrero)

The optimal lambda for variance stabilization is 1.73.

4.12 Time series decomposition

dcmp <- df |>
  model(stl = STL(consumption_gwh)) |>
  components()

dcmp

4.13 Trend and seasonal strength

Feature	Monthly gas consumption in GWh	Monthly mean temperature in Celcius
trend_strength	0.8324117	0.1722841
seasonal_strength_year	0.9463735	0.9527928
seasonal_peak_year	1	7
seasonal_trough_year	8	1
spikiness	51.265714290	0.001106719
linearity	-394.9711522	0.3238554
curvature	-100.9533782	-0.1425342
stl_e_acf1	0.3120111	-0.1448544
stl_e_acf10	0.6141027	0.2663283

References

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). https://OTexts.com/fpp3; OTexts: Melbourne, Australia.