# A Heat-Jarrow-Morton framework for energy markets: a pragmatic approach

Matteo Gardini Department of Mathematics, University of Genoa, Via Dodecaneso
45, 16146, Genoa, Italy, email<EMAIL_ADDRESS>Edoardo Santilli Eni
Plenitude, Via Ripamonti 85, 20136, Milan, Italy, email
<EMAIL_ADDRESS>
###### Abstract

In this article we discuss the application of the Heat-Jarrow-Morton framework
Heath et al. [26] to energy markets. The goal of the article is to give a
detailed overview of the topic, focusing on practical aspects rather than on
theory, which has been widely studied in literature. This work aims to be a
guide for practitioners and for all those who deal with the practical issues
of this approach to energy market. In particular, we focus on the markets’
structure, model calibration by dimension reduction with Principal Component
Analysis (PCA), Monte Carlo simulations and derivatives pricing. As
application, we focus on European power and gas markets: we calibrate the
model on historical futures quotations, we perform futures and spot
simulations and we analyze the results.

Keywords: Stochastic processes, Heat-Jarrow-Morton, Energy Markets, Monte
Carlo, Principal Components Analysis, Calibration, Pricing.

##  1 Introduction

Electricity markets across the world differ for many factors, both
fundamental, such as demand and generation mix, and regulatory. In the most
modern countries one goal in deregulating energy markets is to allow them to
respond to supply and demand variation in a more efficient way. Focusing on US
electricity market, Park et al. [39] has shown that as a result of
deregulation, more competitive and interrelated environments are developing in
the electricity and natural gas markets. In a previous paper Emery and Liu
[22] have discovered that the daily settlement prices of New York Mercantile
Exchange’s (NYMEX’s) California–Oregon Border (COB) and Palo Verde (PV)
electricity futures contracts are cointegrated with the prices of its natural-
gas futures contracts. Such a result has been confirmed by Mjelde and Bessler
[36] which have shown how electricity prices mainly response to shocks in coal
market, whereas Bachmeier and Griffin [3] deeply investigated the level of
market integration between crude oil, coal, and natural gas prices. Moving to
European countries several authors discussed how electricity prices react to
variations in fuel prices. Testing for market integration between natural gas,
electricity and oil prices in the UK in the period in which the natural gas
market was deregulated but not yet linked to the continental European gas
market, Asche et al. [2] have highlighted evidences for an integrated energy
market. Panagiotidis and Rutledge [37] analyzed the relationship between UK
wholesale gas prices and Brent oil price finding co-integration over the
entire sample period(1996–2003). Likewise, using daily price data for Brent
crude oil, NBP UK natural gas and EEX electricity Bencivenga et al. [7] have
shown that gas, oil and electricity markets are integrated. On the other hand
as a result of a robust multivariate long-run dynamic analysis Bosco et al.
[13] have revealed the presence of four highly integrated central European
electricity markets (France, Germany, the Netherlands and Austria). The trend
shared by these four electricity markets appears to be common also to gas
prices, but not to oil prices.

The recent invasion of Ukraine by Russia, and the fear of a possible shortage
in gas supply for Europe, led to an increase in gas and electricity
prices111Such an increase has begun before the invasion, after the Covid-19
pandemic and many complex factors contributed to it. Nevertheless, the effect
of the war on European energy commodities prices has been evident. which has
never seen before as shown in Figure 1. This is easy to explain from an
economical point of view, since in Europe natural gas is used to produce the
$19.2\%$ of electricity and natural gas power plants usually play the role of
marginal technology in the electricity supply curve.

(a) Future prices, delivery year 2023

(b) Spot prices

Figure 1: German power and natural TTF gas daily prices from1/1/2020 to the
31/12/2021.

At this point it should be evident that integration between energy markets
must be taken into account if one is interested in energy commodities
modeling, risk management or in derivatives pricing.

Over the years many approaches have been proposed to model energy markets in a
univariate setting. Pioneering papers in this field dates back at the
beginning of the century: Schwartz and Smith [47], Lucia and Schwartz [32] and
Schwartz [46] focused mainly purely Gaussian framework, whereas Cartea and
Figueroa [19] proposed a mean-reverting model with jumps and a deterministic
seasonality component for electricity spot prices. Saifert and Uhrig-Homburg
[42] compared different modeling approaches in power markets, whereas a good
summary of energy markets modeling is contained in Benth et al. [9]. Many non-
Gaussian models based on Lévy processes which have been proposed for equities,
such as the variance gamma (Madan et al. [34]), the jump-diffusion model
(Merton [35]) and the normal-inverse Gaussian process (Barndorff-Nielsen [5]):
in particular a Two-Factor version of this process have been recently applied
to the energy context by Piccirilli et al. [41]. Furthermore, many stochastic
volatility models, like the ones proposed by Heston [27] and Bates [6], can be
adapted to model commodity prices behavior. All these powerful tools can be
used to properly consider many stylized facts such as jumps in price
trajectories, skewness and fat-tails in log-returns distribution and
volatility smiles. A review of financial modeling with jump processes can be
found in Cont and Tankov [20].

Financial modeling in univariate setting has been deeply investigated, but
challenging issues arise when we scale to a multi-commodities market. Within
this context, the former modeling techniques become harder to apply in
practice and literature is not as rich as in the one dimensional framework.
Petroni and Sabino [40] have shown how some standard models such as the ones
proposed by Black and Scholes [11], Schwartz and Smith [47] and Cartea and
Figueroa [19] can be extended to a multivariate context by adding dependent
jumps which are modeled using self-decomposable laws, whereas Kiesel and
Kusterman [30] introduced a structural model to properly consider the market
coupling effect in electricity spot markets, in the spirit of what has been
proposed by Carmona and Coulon [18].

A widely recognized approach for energy markets modeling has been proposed by
Benth and Saltyte-Benth [8] which adapted the framework introduced by Heath et
al. [26] to energy futures market. This modeling technique, together with the
calibration of the underlying model, has been studied by many authors such as
Sclavounos and Ellefsen [48], Hinderks et al. [28], Benth et al. [10],
Broszkiewicz-Suwaj and Aleksander [15], Edoli et al. [21] and Feron and Gruet
[23]. Despite their mathematical completeness and accurateness, these articles
seems hard to be used in practice since many practicalities are not covered.
The data preparation, a clear explanation on how the model should be
implemented in practice, a “practitioner” interpretation of the results, the
management of typical issues arising during the implementation are often
missing. This works aims at filling this gap, by collecting all the results
presented so far concerning the application of the Heat-Jarrow-Morton (HJM)
framework to energy markets. By focusing on European power and gas futures
markets, we present a very general approach for data analysis and preparation,
for model calibration and, finally, for simulation. Furthermore, we discuss
some model limitations and we suggest some possible extensions, trying to
preserve both the numerical and the mathematical tractability of the
framework.

The article is organized as follows: Section 2, focuses on a slice of European
power and gas futures markets, and shows how the HJM approach might be the
appropriate modeling framework. Section 3 introduces the model and briefly
explain its behavior with the support of some “toy-examples”. In this section
we also show how the PCA can be used in order to calibrate the market
parameters. Section 4 considers real-market data, deals with different
approaches to data preparation and shows how to properly calibrate the model
on power and gas European futures prices. In Section 5 we show simulation
results of the whole market considered in Section 4: we analyze the outcomes
and we briefly discuss model’s strengths and limitations. Section 6 concludes
the work and discusses some possible model extensions.

##  2 Market structure and analysis

In this section we focus on the European power and gas futures markets. The
results we obtain are valid within this framework, but the approach can be
easily adapted to any market. In particular, we consider data from the
European Energy Exchange (EEX) for the power markets, whereas data regarding
the natural gas markets comes from the Intercontinental Exchange (ICE).

Recently, the power and gas futures markets in Europe has experimented an era
of expansion and increasing interest. Within this market many contracts are
present but the most traded ones are those with monthly, quarterly and yearly
delivery. For example, a power future 2024 calendar is a contract between two
counterparts to buy or sell a specific volume of energy in MWh at fixed price,
decided at time $t$, for all the hours of the year. In this case
$T^{s}=1/1/2024$ and $T^{e}=31/12/2024$. Moreover, at time $T^{s}$ the product
expires and hence $T^{s}$ plays the role of maturity222Actually the product
with delivery starting at $T^{s}$ is traded until the day before the beginning
of the delivery period, but in order to simplify the notation we consider
$T^{s}$ as the maturity.. Futures contracts with monthly or quarterly delivery
are defined similarly. We will denote the price of such a contract at time $t$
by $F(t,T^{s},T^{e})$. Concerning the gas market, the situation is slightly
different, but for the aim of this paper we can consider it to be similar to
power one. In particular, we can assume that the same contract
$F(t,T^{s},T^{e})$ written on gas, will delivery the gas you need to produce 1
MWh of electricity for the whole delivery period $\left[T^{s},T^{e}\right]$.
On the other hand, we denote by $F(t,T)$, the so called _fixed delivery_
futures contract: $F(t,T)$ is the price at time $t$ of a contract with
delivery $T$, for example one year,one month, three months and so on from now.
These latter products are not quoted in the market but can be easily obtaining
by rearranging $F(t,T^{s},T^{e})$ as we will show in Section 4.1. From a
modeling point of view great simplification arises by defining the dynamic of
$F(t,T)$ and hence obtaining that of $F(t,T^{s},T^{e})$.

Without loss of generality, in this paper we focus on four power markets and
on two natural gas markets. In particular we consider the German (DE), Italian
(IT), French (F7) and Swiss (CH) power futures markets and the TTF and PSV
which are the Dutch and the Italian hubs for the natural gas, respectively. Of
course, the analysis can be extended to an arbitrary number of markets.
The final goal is to model the dynamic of the whole forward curve for the
aforementioned markets under a HJM framework. We can assume that each of
traded contract acts as a source of uncertainty (which we call a “random
factor”) for the determination of the forward curve dynamic. On the other
hand, as we stated before, all these markets are co-integrated and hence the
hope is that we can use a small number of random factors to successfully model
the forward curves dynamic. In oder to verify this assumption, we consider
daily futures prices between the $1^{st}$ of January 2020 to the $31^{st}$ of
December 2022 and we compute the correlation of the log-returns, following the
approach proposed by Sclavounos and Ellefsen [48]. In Figure 2 we plot the
correlation surfaces between daily log-returns of several futures products
with fixed delivery333As mentioned before, in energy markets it is a common
practice to distinguish between general future contracts and those ones with
fixed delivery. Fixed delivery products refer to delivery occurring after a
given amount of time with respect to the present date. For example, if today
is the 19th of November 2023, $F(t,T_{1})$ with $T_{1}=\frac{1}{12}$ refers to
the to product with delivery on the next month, namely December 2023: it is
customary to refer to it as $M1$ (a short hand for month plus one). An example
will make the distinction clearer. $F(t,T^{s},T^{e})$ refers to an absolute
delivery, for example December 2023 and hence, $T^{s}=1/12/2023$ and
$T^{e}=31/12/2023$. If today is the 19th October 2023, this product will be
the $M2$ and will be denoted by $F(t,T_{2})$, with $T_{2}=\frac{2}{12}$. As
time goes on if we are on the 19th November 2023 we observe the same product
$F(t,T^{s},T^{e})$ but this time this product correspond to the $M1$ since its
delivery is only one month far from now. We denote this product by
$F(t,T_{1})$ with $T_{1}=\frac{1}{12}$.: we observe that the linear
correlation coefficient is significant and hence this bodes well that only few
stochastic factors drive the whole structure of the forward curve.

Figure 2: Log-returns correlation surfaces for different commodities.

Most financial models relies upon stochastic processes in continuous time and,
mainly, on Lévy processes: the well known Brownian motion is just the simplest
of them. Working in a Lévy framework leads to the development models
characterized by a very reach structure which can be used to efficiently
include many stylized facts. The interested reader can refer to Sato [44] and
Applebaum [1] for an overview on Lévy processes and to Cont and Tankov [20]
for applications to financial markets. One of the assumptions of Lévy
processes is the independence of increments. For this reason, before using
Lévy processes for modeling purposes, one has to check that the increments are
independent. Following the approach proposed in Brigo et al. [14], we compute
the auto-correlation function (ACF) on six different time series, one for each
market, on the calendar product with delivery the year 2023. We consider the
series of daily log-returns $x_{1},x_{2},\dots,x_{n}$, with

$x_{i}=\ln\frac{S_{t_{i+1}}}{S_{t_{i}}},$

where $S_{t_{i}}$ denotes the price of a given risky asset at time $t_{i}$ and
we compute the ACF with lag $k$ as:

$ACF(k)=\frac{1}{(n-k)\hat{v}}\sum_{i=1}^{n-k}\left(x_{i}-\hat{m}\right)\left(x_{i+k}-\hat{m}\right),\quad
k=1,\dots,20,$

where $\hat{m}$ and $\hat{v}$ are the sample mean and variance. Roughly
speaking we can consider the ACF as an estimate of the correlation between the
random variables $X(t_{i})$ and $X(t_{i+k})$. ACF for the six selected
products above are shown in the charts in Figure 3. For the all of them we do
not observe any significant lags in the historical return time series, which
means the independence assumption is acceptable in this case. Changing the
delivery period of the product we get similar results. Therefore Lévy
processes can be used in order to properly model the futures prices.

Figure 3: Sample ACF computed on log-returns.

Within the classical HJM framework log-returns are assumed to be normally
distributed. As observed by many authors in most financial markets log-returns
are not normally distributed, but their distribution is often skewed and
presents heavy tails or fat tails effect. Furthermore, the log-returns
volatility is not constant but often clusters appear. In Figure 4 we plot the
empirical probability density function of log-returns compared with the normal
one fitted on the same set of data as above. In all cases we observe that the
real log-returns distribution is peaked and presents tails which are heavier
than the normal distribution ones. Furthermore, all distributions result to be
skewed.

In light of this results, by choosing a normal distribution we might lose some
market peculiarities. On the other hand, in order to get a simple and stable
calibration methodology the hypothesis of normality in log-returns is commonly
accepted by practitioners. Gaussianity hypothesis in log-returns might be
relaxed by using Lévy processes or Lévy copulas, as proposed by Panov and
Samarin [38] and Cont and Tankov [20], but this approach is very hard to
handle from a practical point of view. Indeed, the calibration step is hard to
tackle, especially if the number of underlying asset is large. Several
authors, such as Luciano and Semeraro [33], Schoutens [45], Ballotta and
Bonfiglioli [4], Buchmann et al. [16] and Gardini et al. [24], investigated
other techniques in order to consider dependence in log-returns remaining in a
Lévy framework. Most of them work quite well if the number of the risky assets
to remains small, but complications arise when one deals with many risky
underlying assets since the number of parameters rapidly grows as the number
of underlying assets increases. Consequently calibration becomes hard to
perform in practice and its results might be unreliable. For these reasons,
the Gaussian framework remains a milestone among practitioners in multi
commodity energy markets. On the other hand, if a focus on a single product is
required, a more general model among the ones we listed should be considered.

Figure 4: Daily futures log-returns densities for calendar 23 products.

In the next sections we apply the HJM framework to energy markets, with a
particular focus on European power and natural gas markets. Nevertheless, the
approach is very general and could be easily adapted to other commodity
markets such as oil, precious metals and agricultural products.

##  3 The model

In this section we discuss in detail the HJM framework applied to energy
markets. In particular we start from some “toy-models” which are useful to fix
the main modeling concepts. Once that the dynamics and the calibration
procedures for these simpler models are clear, we introduce the most general
framework which turns out to be a simple extension of the previous setting.

###  3.1 A single factor toy-model

Consider the price of the fixed delivery future contract, signed at $t_{0}$,
$F(t,T)$ at time $t\in\left[t_{0},T\right]$ for a fixed delivery $T$. Assume
that we have only a single source of uncertainty and, assuming to work under
the risk neutral measure $\mathbb{Q}$, consider a dynamic of the following
type:

$\frac{dF(t,T)}{F(t,T)}=\sigma(t,T)dW(t),\quad F(t_{0},T)=F(0,T)\;a.s.,$ (1)

where the volatility $\sigma(t,T)$ is assumed to be a deterministic time
dependent function $\sigma(t,T):\left[t_{0},T\right]\mapsto\mathbb{R}^{+}$ and
$F(0,T)\in\mathbb{R}^{+}$ is the value of the future contract at time $t_{0}$.

According to Samuelson [43] the term structure of commodity forward price
volatility typically declines with contract horizon: this is what is commonly
known as Samuelson’s effect. Therefore, it is customary to assume a volatility
function which depends on the time to maturity $T-t$, namely
$\sigma\left(t,T\right)=\sigma(T-t)$. In particular, $\sigma(T-t)$ will be
decreasing in $T-t$, reflecting the fact that contract with longer maturities
are less volatile than contracts with a shorter one.

Using very basic Itō’s calculus we can solve the Equation (1), obtaining

$F(t,T)=F(t_{0},T)\exp\left\\{-\frac{1}{2}\int_{t_{0}}^{t}\sigma\left(T-s\right)^{2}ds+\int_{t_{0}}^{t}\sigma(T-s)dW(s)\right\\}.$

The existence of an explicit solution for the stochastic differential Equation
(1) is extremely important in order to simulate the trajectories of process
$F(T)=\left\\{F(t,T);t_{0}\leq t\leq T\right\\}$ exactly. Some possible
realizations of the process over one year time horizon for a volatility
function of the form:

$\sigma(T-t)=0.8e^{-2(T-t)},$

are shown in Figure 6. Observe that, at the beginning the volatility is low,
since the time to maturity is large, whereas as $T-t\to 0$ the volatility
increases, according to the Samuelson’s effect.

Figure 5: Annualized volatility of the toy model including the Samuelson’s
effect. Figure 6: Possible simulations of the forward prices $F(t,T)$.

When one performs simulations, “sanity checks” are important. By “sanity
check” we mean the comparison of a numerical computation with a theoretical
one. For example, by Itō’s isometry, it is very easy to check that:

$Var\left[\ln
F(t,T)\right]=Var\left[\int_{t_{0}}^{t}\sigma(T-s)dW(s)\right]=\int_{t_{0}}^{t}\sigma\left(T-s\right)^{2}ds,\;t\in\left[t_{0},T\right].$
(2)

Hence, one can estimate numerically the variance of $F(t,T)$ for each
$t\in\left[t_{0},T\right]$ and compare this result with the one given by
Equation (2). From Figure 7 we observe that the numerical quantities and the
theoretical ones are very close and hence we are guaranteed that the numerical
simulation scheme is correctly implemented.

Figure 7: Sanity check of the volatility of log-prices in the toy-model by
using $5000$ simulations.

The knowledge of and explicit expression for variance as the one in Equation
(2) is important also for option pricing. Following Böerger et al. [17], if we
assume the following dynamic for the forward price:

$\frac{dF(t,T)}{F(t,T)}=\boldsymbol{\sigma}(t,T)\cdot d\boldsymbol{W}(t),$

where
$\boldsymbol{\sigma}(t,T)=\left(\sigma_{1}(t,T),\dots,\sigma_{n}(t,T)\right)$
and $\boldsymbol{W}=\left(W_{1},\dots,W_{n}\right)$ is a $n$-dimensional
standard Brownian motion with independent components, the price of a call
option with maturity $T_{0}\leq T$ and strike price $K$ is given by the
standard Black formula:

$C(t_{0},T_{0},K)=e^{-r(T_{0}-t_{0})}\left(F(t_{0},T)\mathcal{N}\left(d_{1}\right)-K\mathcal{N}\left(d_{2}\right)\right),$
(3)

where $r\geq 0$ is the risk-free rate, $\mathcal{N}\left(\cdot\right)$ is the
cumulative distribution function of a standard normal random variable and:

$\displaystyle d_{1}$
$\displaystyle=\frac{\ln\frac{F(T_{0},T)}{K}+\frac{1}{2}Var\left[\log
F(T_{0},T)\right]}{\sqrt{Var\left[\log F(T_{0},T)\right]}}$ $\displaystyle
d_{2}$ $\displaystyle=d_{1}-\sqrt{Var\left[\log F(T_{0},T)\right]},$

where $Var\left[\log F(T_{0},T)\right]$ can be easily computed by Itō’s
isometry and depends on the form we chose for $\boldsymbol{\sigma}(T-t)$.
For example, if we consider a bi-dimensional standard Brownian motion and we
assume that:

$\displaystyle\sigma_{1}(t,T)$ $\displaystyle=0.8e^{-2(T-t)},$
$\displaystyle\sigma_{2}(t,T)$ $\displaystyle=0.2,$

we obtain:

$Var\left[\log
F(T_{0},T)\right]=\frac{0.64}{4}\left[e^{-4(T-T_{0})}-e^{-4(T-t_{0})}\right]+0.04\left(T_{0}-t_{0}\right).$

The price of European call option with different strike prices $K$ are shown
in Figure 8

Figure 8: Pricing of European call options with different strike prices. Exact
formula is compared to a Monte Carlo pricing with $10^{7}$ simulations.

Of course, any form can be assumed for $\sigma(T-t)$, but, from a practical
point of view a step-wise volatility structure is usually considered. Indeed,
if we consider monthly fixed delivery, such as $M0,M1,\dots,M36$, with $M=36$
delivery dates $\left\\{T_{j}\right\\}_{j=1}^{M}$, the volatility function
assumes the following form:

$\sigma(T-t)=\sum_{j=1}^{M}\sigma_{j}\mathbbm{1}_{I_{j}}(T-t),$

where $\sigma_{j}\in\mathbb{R}^{+}$ and $I_{j}=\left(T_{j-1},T_{j}\right]$ and
where $T_{j}$ is the delivery of the $j$-th contract. As mentioned before, we
typically do not observe monthly products but quarters, season and yearly
contract. On the other hand it is possible to extract from the available
market products the monthly quotations. The data preparation methodology and
the calibration algorithm we use to estimate $\sigma_{j}$ is discussed in
Section 4.

###  3.2 A model for a multi-commodity market

Now we extend the “toy-model” we presented in the previous section to a multi-
commodity framework. Assume we have $k=1,\dots,K$ markets (DE, F7, IT, TTF and
so on), each of them consist in the same number $M_{k}=M$ for $k=1,\dots,K$ of
fixed delivery futures contracts $F^{k}(t,T_{m})$. Recall that a contract
which delivers at the beginning of the next month, independently of the
trading date is stated by $M1$, which means that the delivery is the next
month after the current date. As mentioned above, these products are not
quoted on the market, where only futures contracts $F(t,T^{s},T^{e})$ can be
observed, but they can be easily computed from the quoted ones. As a first
step of our modeling framework we start by defining the dynamic of the fixed
delivery futures contract. Assume that each of the the contracts is a “random
factor” which might have an impact on the whole market dynamic. Hence we have
a total number of $\tilde{N}=M\cdot K$ random factors. The dynamic of a future
product $k=1$ with fixed delivery $T_{m}$ is given by:

$\frac{dF^{k}(t,T_{m})}{F^{k}(t,T_{m})}=\sum_{j=1}^{\tilde{N}}\sigma_{kj}dW_{j}(t),\quad
F^{k}(t_{0},T_{m})=F^{k}(0,T_{m})\;a.s.,\;t\geq 0.$ (4)

where $W_{j}=\left\\{W_{j}(t);t\geq 0\right\\}\;j=1,\dots,\tilde{N}$ are
independent Brownian motions. From Equation (4) observe that the dynamic of
the single monthly product $F^{k}(t,T_{m})$ potentially depends on the ones of
all other monthly futures products. From empirical evidence, energy markets
are strongly co-integrated and hence it is clear that considering $\tilde{N}$
independent Brownian motions appears to be unreasonable. A possible
methodology, based on the dimensional reduction inherited from the Principal
Component Analysis (PCA), to select a lower number of stochastic factors, will
be presented in Section 4.
Observe that Equation (4) can be written in a matrix form as:

$\frac{d\boldsymbol{F}(t,T)}{\boldsymbol{F}(t,T)}=\boldsymbol{\sigma}\cdot
d\boldsymbol{W}(t),$

where $\boldsymbol{W}=\left(W_{1},\dots,W_{\tilde{N}}\right)$ is a standard
Brownian motion with independente components, $\boldsymbol{\sigma}$ is a
$\tilde{N}\times\tilde{N}$ matrix of the form:

$\boldsymbol{\sigma}=\begin{matrix}\boldsymbol{\sigma}^{1}\left\\{\vphantom{\begin{matrix}\sigma_{1,1}\\\
\sigma_{2,1}\\\ \dots\\\ \sigma{M,1}\\\ \hline\cr\end{matrix}}\right.\\\
\boldsymbol{\sigma}^{2}\left\\{\vphantom{\begin{matrix}\sigma_{M+1,1}\\\
\sigma_{M+2,1}\\\ \dots\\\ \sigma{2M,1}\\\ \hline\cr\end{matrix}}\right.\\\
\boldsymbol{\sigma}^{k}\left\\{\vphantom{\begin{matrix}\vdots\\\ \vdots\\\
\vdots\\\ \hline\cr\end{matrix}}\right.\\\
\boldsymbol{\sigma}^{K}\left\\{\vphantom{\begin{matrix}\sigma_{M(K-1)+1,1}\\\
\sigma_{M(K-1)+2,1}\\\ \dots\\\ \sigma{MK,1}\\\ \end{matrix}}\right.\\\
\end{matrix}\begin{bmatrix}\sigma_{1,1}&\sigma_{1,2}&\sigma_{1,3}&\dots&\sigma_{1\tilde{N}}\\\
\sigma_{2,1}&\sigma_{2,2}&\sigma_{M+2,3}&\dots&\sigma_{2,\tilde{N}}\\\
\dots&\dots&\dots&\dots&\dots\\\
\sigma_{M,1}&\sigma_{M,2}&\sigma_{M,3}&\dots&\sigma_{2M,\tilde{N}}\\\
\hline\cr\sigma_{M+1,1}&\sigma_{M+1,2}&\sigma_{M+1,3}&\dots&\sigma_{M+1,\tilde{N}}\\\
\sigma_{M+2,1}&\sigma_{M+2,2}&\sigma_{M+2,3}&\dots&\sigma_{M+2,\tilde{N}}\\\
\dots&\dots&\dots&\dots&\dots\\\
\sigma_{2M,1}&\sigma_{2M,2}&\sigma_{2M,3}&\dots&\sigma_{2M,\tilde{N}}\\\
\hline\cr\vdots&\vdots&\vdots&\vdots&\vdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
\hline\cr\sigma_{M(K-1)+1,1}&\sigma_{M(K-1)+1,2}&\sigma_{M(K-1)+1,3}&\dots&\sigma_{M(K-1)+1,\tilde{N}}\\\
\sigma_{M(K-1)+2,1}&\sigma_{M(K-1)+2,2}&\sigma_{M(K-1)+2,3}&\dots&\sigma_{M(K-1)+2,\tilde{N}}\\\
\dots&\dots&\dots&\dots&\dots\\\
\sigma_{MK,1}&\sigma_{MK,2}&\sigma_{MK,3}&\dots&\sigma_{MK,\tilde{N}}\end{bmatrix}$

where $\boldsymbol{\sigma}^{k}\in\mathbb{R}^{M\times\tilde{N}}$ denotes the
matrix associated to random factors of $k$-th market product.

Once we have defined the dynamic for the fixed delivery products that of the
$F(t,T^{s},T^{e})$ is a direct consequence. Assume now that $T^{s}$ is an
absolute date, for example $31/12/2023$. Then, the dynamic of
$F^{k}(t,T^{s},T^{e})$ for $t\in\left[t_{0},T^{s}\right]$ is given by
444Observe that we are assuming that $\forall T^{s}$ there exists
$m_{0}\in\left[1,\dots,M\right]$ such that $T_{s}=T_{m_{0}}$. For example
today $t_{0}=19/2/2023$ we observe $F(t,T^{s},T^{e})$ with $T^{s}=1/5/2023$
and $T^{e}=31/5/2023$ we are assuming modeling the product
$F(t_{0},T_{m_{0}})$ in such a way that $T^{s}-t_{0}=T_{m_{0}}-t_{0}$.:

$\frac{dF^{k}(t,T^{s},T^{e})}{F^{k}(t,T^{s},T^{e})}=\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\hat{\sigma}_{ij}^{k}\mathbbm{1}_{I_{i}}(T-t)dW_{j}(t).$
(5)

The explicit solution of the Equation (5) for $t\in\left[t_{0},T^{s}\right]$
is given by:

$\begin{split}F^{k}(t,T^{s},T^{e})=F^{k}(t_{0},T^{s},T^{e})\exp\left\\{-\frac{1}{2}\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\left(\sigma_{ij}^{k}\right)^{2}\int_{t_{0}}^{t}\mathbbm{1}_{I_{i}}\left(T^{s}-s\right)ds\right.\\\
\left.+\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\sigma_{ij}^{k}\int_{t_{0}}^{t}\mathbbm{1}_{I_{i}}\left(T^{s}-s\right)dW_{j}(s)\right\\}.\end{split}$
(6)

The existence of a close form solution allows us to simulate the paths of the
process in an exact way, without using any discretization method, such as the
Euler’s of Millstein’s one (see Seydel [49]). Once again, observe that if
$t\approx t_{0}$ then $T^{s}-t$ is “large” and hence we are considering the
volatility of the products with a large maturity. As soon as time $t$ goes on,
the time to delivery $T^{s}-t$ gets smaller and smaller hence we are taking
into account the volatility of the fixed delivery products with a short time
to maturity.

Also in this case a sanity check can be performed. Assume we simulate the
process following Equation (6) for a very short period, namely $T^{s}-t\approx
T^{s}-t_{0}$. Hence Equation (6) simplifies to:

$F^{k}(t,T_{m_{0}})=F^{k}(t_{0},T_{m_{0}})\exp\left\\{-\frac{1}{2}\Delta
t\sum_{j=1}^{\tilde{N}}\left(\sigma_{m_{0}j}^{k}\right)^{2}+\sum_{j=1}^{\tilde{N}}\sigma_{m_{0}j}^{k}W_{j}(t)\right\\},$
(7)

where $\Delta t=t-t_{0}$ and $m_{0}$ is such that $T_{m_{0}}=T^{s}$. By
computing the log-return $x^{k}(t)=\ln F(t,T_{m_{0}})-\ln F(t_{0},T_{m_{0}})$
we get:

$x^{k}(t)=-\frac{1}{2}\Delta
t\sum_{j=1}^{\tilde{N}}\left(\sigma_{m_{0}j}^{k}\right)^{2}+\sum_{j=1}^{\tilde{N}}\sigma_{m_{0}j}^{k}dW_{j}(t),$

and if we compute its variance we get:

$Var\left[x^{k}(t)\right]=\Delta
t\sum_{j=1}^{\tilde{N}}\left(\sigma_{m_{0}j}^{k}\right)^{2}.$

Hence, we have a simple way to check that the variance of the simulated
product is correct, by summing up the squares of the entries of the matrix
$\boldsymbol{\sigma}^{k}$ and multiplying by $\Delta t$.
Clearly if the approximation $T_{m_{0}}-t\approx T-t_{0_{0}}$ does not hold,
namely for larger $t$ such expression is no longer valid. On the other hand,
simple computations show that:

$Var\left[\log
F^{k}(t,T^{s},T^{e})\right]=\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\left(\sigma_{ij}^{k}\right)^{2}\int_{t_{0}}^{t}\mathbbm{1}_{I_{i}}\left(T^{s}-s\right)ds,$
(8)

which can be easily used to check the correctness of the simulations.

As stated in Section 3.1 once that the expression of $Var\left[\log
F^{k}(t,T^{s},T^{e})\right]$ is known, a Black’s style pricing formula of the
form shown for European call options as the one shown in Equation 3 can be
easily derived.

###  3.3 The spot dynamic

Starting from Equation (6) we can retrieve the dynamic for the spot prices by
defining the spot prices $S(t)$ as:

$S(t)=\lim_{\begin{subarray}{c}T^{s}\to t\\\ T^{e}\to
t\end{subarray}}F(t,T^{s},T^{e}).$

Passing to the limit in Equation (6) we have that:

$\begin{split}S^{k}(t)=F^{k}(t_{0},t)\exp\left\\{-\frac{1}{2}\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\left(\sigma_{ij}^{k}\right)^{2}\int_{t_{0}}^{t}\mathbbm{1}_{I_{i}}\left(t-s\right)ds\right.\\\
\left.+\sum_{j=1}^{\tilde{N}}\sum_{i=1}^{M}\sigma_{ij}^{k}\int_{t_{0}}^{t}\mathbbm{1}_{I_{i}}\left(t-s\right)dW_{j}(s)\right\\}.\end{split}$
(9)

Observe that it is important to understand how volatility behaves as time $t$
increases. If $t\approx t_{0}$ then the volatility of the spot price
$S^{k}(t)$ the one of the fixed delivery futures products with the shorter
time to maturity. As time $t$ increases, the spot price $S^{k}(t)$ somehow
includes all the volatility effects from the products with shorter time to
maturity to the ones with a longer one. This is reasonable from an economical
point of view. Indeed, if we want to simulate the process $S^{k}(t)$ in a year
it would have been the $M12$, the $M11$ and so on up to $M0$ and hence it
includes all their volatility. Moreover, an expression for
$Var\left[S^{k}(t)\right]$ similar to the one in Equation (8) can be easily
obtained.

From a practical point of view, especially for power prices, it is customary
to consider a hourly granularity. $F^{k}(t_{0},t)$ represents the power
forward price today for hour $t$ and hence the function $F^{k}(t_{0},t)$ for
$t\in\left[t_{0},T\right]$ for a fixed $T$, represents the hourly forward
curve for the product $k$. In simulations routines, in order to simplify the
simulation, we assume the the “shocks” with respect to the hourly forward
curve are daily. For commodities with daily granularity, such as the natural
gas, everything is performed on a daily basis. The hourly or daily forward
curve should be inferred from futures market and can be obtained by using
different methodologies, such as the one proposed by Benth et al. [9, Chapter
7].

###  3.4 A multidimensional toy example

In order to make things clearer, in this section we present a toy version of
the model we presented above. Assume we have only four products with different
deliveries $T_{m},\;m=1,\dots,M$ with $M=4$. Assume we are focusing on a
single market (for example the power futures market in Germany) and hence
$K=1$. The fixed delivery products have the following dynamic:

$\frac{dF\left(t,T_{m}\right)}{F(t,T_{m})}=\sigma_{m1}dW_{1}(t)+\sigma_{m2}dW_{2}(t)+\sigma_{m3}dW_{3}(t)+\sigma_{m4}dW_{4}(t)=\sum_{j=1}^{\tilde{N}}\sigma_{mj}dW_{j}(t).$

Observe that the components of the multidimensional standard Brownian motion
$\boldsymbol{W}=\left\\{\left(W_{1}(t),W_{2}(t),W_{3}(t),W_{4}(t)\right);t\in\left[t_{0},T\right]\right\\}$,
can be both correlated and independent. The dynamic in Equation (3.4) can be
written in a matrix form as:

$\frac{d\boldsymbol{F}(t,T)}{\boldsymbol{F}\left(t,T\right)}=\boldsymbol{\sigma}\cdot
d\boldsymbol{W}(t),\quad t\in\left[t_{0},T\right]$

where :

$\begin{split}\boldsymbol{F}(t,T)=\begin{bmatrix}F\left(t,T_{1}\right)\\\
F\left(t,T_{2}\right)\\\ F\left(t,T_{3}\right)\\\
F\left(t,T_{4}\right)\end{bmatrix},\qquad
T=\min_{m\in\left[1,M\right]}T_{m},\qquad\boldsymbol{\sigma}=\begin{bmatrix}\sigma_{11}&\sigma_{12}&\sigma_{13}&\sigma_{14}\\\
\sigma_{21}&\sigma_{22}&\sigma_{23}&\sigma_{24}\\\
\sigma_{31}&\sigma_{32}&\sigma_{33}&\sigma_{34}\\\
\sigma_{41}&\sigma_{42}&\sigma_{43}&\sigma_{44}\\\ \end{bmatrix}\end{split}$

In the market you do not directly observe the matrix $\boldsymbol{\sigma}$ but
it can be efficiently estimated from real market data. In order to do so, we
assume that $N$ equally spaced daily observation at times $t_{0}\leq
t_{1}\leq\dots\leq t_{N}$ of the process $\boldsymbol{F}$ are given and let
$\Delta t=t_{i+1}-t_{i}$ defined in fraction of years, for example $\Delta
t=\frac{1}{260}$. Define the log-return for the product with delivery $T_{m}$
as:

$X_{i}^{m}=\ln\frac{F(t_{i+1},T_{m})}{F\left(t_{i},T_{m}\right)}$

and assume that the vector
$\boldsymbol{X}=\left[X^{1},X^{2},X^{3},X^{4}\right]$ is normally distributed
with mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$. Since
$\boldsymbol{\Sigma}$ is a covariance matrix it is symmetric and positive-
definite matrix and hence it factorizes as:

$\boldsymbol{\Sigma}=\boldsymbol{C}\boldsymbol{\Gamma}\boldsymbol{C}^{T}=\boldsymbol{C}\boldsymbol{\Gamma}^{\frac{1}{2}}\left(\boldsymbol{C}\boldsymbol{\Gamma}^{\frac{1}{2}}\right)^{T}$

for $\boldsymbol{\Gamma}\in\mathbb{R}^{M\times M}$ diagonal matrix and
$\boldsymbol{C}\in\mathbb{R}^{M\times M}$ is an orthogonal matrix such that
$\boldsymbol{C}^{T}\boldsymbol{C}=\boldsymbol{I}$, where $\boldsymbol{I}$ is
the identity matrix. On the other hand we have that
$\boldsymbol{\Sigma}=\boldsymbol{X}^{T}\boldsymbol{X}$ and it can be proved
that:

$\boldsymbol{\Sigma}=\Delta t\boldsymbol{\sigma}\boldsymbol{\sigma}^{T}.$

Therefore, $\boldsymbol{\sigma}$ can be estimated as:

$\boldsymbol{\sigma}=\frac{\boldsymbol{C}\boldsymbol{\Gamma}^{\frac{1}{2}}}{\Delta
t^{\frac{1}{2}}}.$

For this reason we can easily estimate the matrix $\boldsymbol{\sigma}$ from
the observed data using Algorithm 1.

Algorithm 1 Estimation of $\boldsymbol{\sigma}$

1:Assume to observe $N$ realization of the random vector $\boldsymbol{X}$ and
list them into the matrix $\hat{\boldsymbol{X}}$.

2:Compute the sample covariance matrix $\hat{\boldsymbol{\Sigma}}$ as
$\hat{\boldsymbol{\Sigma}}=\hat{\boldsymbol{X}}^{T}\hat{\boldsymbol{X}}$.

3:Diagonalize the matrix
$\hat{\boldsymbol{\Sigma}}=\boldsymbol{C}\boldsymbol{\Gamma}\boldsymbol{C}^{T}$.

4:$\hat{\boldsymbol{\sigma}}$ is given by:
$\hat{\boldsymbol{\sigma}}=\boldsymbol{C}\boldsymbol{\Gamma}^{1/2}\Delta
t^{-1/2}$.

######  Remark 1.

Observe that the matrix $\hat{\boldsymbol{\sigma}}$ defined in Algorithm 1 is
not unique, but it is unique up to unitary transformation. Indeed assume that
$\boldsymbol{U}\in\mathbb{R}^{M\times M}$ is such that
$\boldsymbol{U}\boldsymbol{U}^{T}=\boldsymbol{I}$. Then:

$\boldsymbol{\Sigma}=\Delta t\boldsymbol{\sigma}\boldsymbol{\sigma}=\Delta
t\boldsymbol{\sigma}\boldsymbol{U}\boldsymbol{U}^{T}\boldsymbol{\sigma}^{T}=\Delta
t\boldsymbol{\sigma}\boldsymbol{U}\left(\boldsymbol{\sigma}\boldsymbol{U}\right)^{T}$

and hence we can define
$\bar{\boldsymbol{\sigma}}=\boldsymbol{\sigma}\boldsymbol{U}$ and this is a
$M\times M$ matrix with the property that:

$\boldsymbol{\Sigma}=\Delta
t\bar{\boldsymbol{\sigma}}\bar{\boldsymbol{\sigma}}^{T}.$

Within our toy example assume that:

$\boldsymbol{\sigma}=\begin{bmatrix}0.15&0.019&-0.13&0.018\\\
0.25&0.014&-0.19&0.015\\\ 0.185&0.012&-0.13&0.018\\\
0.125&0.044&-0.131&0.043\\\ \end{bmatrix},$

and hence simulate the trajectories of the process $\boldsymbol{F}$ for a
given time horizon $\left[t_{0},T\right]$. If we plot the trajectories
obtained by simulating from in Equation (1) we obtain the results in Figure 9.
Hence compute the log-returns and get $\hat{\boldsymbol{\sigma}}$ following
Algorithm 1, which gives:

$\hat{\boldsymbol{\sigma}}=\begin{bmatrix}-0.0004&-0.0056&-0.0054&0.1964\\\
0.0003&-0.0037&0.0203&0.3090\\\ -0.0002&0.0087&0.0141&0.2228\\\
0.0001&0.0017&-0.0462&0.1809\\\ \end{bmatrix}.$

Depending on the number of time step you choose and on the seed of random
number generator you could obtain a different matrix
$\hat{\boldsymbol{\sigma}}$.

Figure 9: Possible simulation of the forward prices process $\boldsymbol{F}$ .

From Figure 9 It appears clear that the whole market can be properly described
by using less than four risk factors. We get the same conclusion by looking at
the matrix correlation $\rho$ of log-returns and we obtain:

$\rho=\begin{bmatrix}1.00&0.99&0.99&0.98\\\ 0.99&1.00&0.99&0.95\\\
0.99&0.99&1.00&0.95\\\ 0.98&0.95&0.95&1.00\\\ \end{bmatrix},$ (10)

If we perform the PCA analysis we get that the eigenvalues of the covariance
matrix $\hat{\boldsymbol{\Sigma}}$ are:

$\lambda=\left[0.8325,0.0107,0.0005,0.0000\right]\cdot 10^{-3}.$

It turns out that the first two principal components (eigenvectors) associated
to the first two eigenvalues explains more than the $99\%$ of the variance.
For this reason we can use only two random (i.e. Brownian motion) factor to
properly model all the market, instead of the original four. Define:

$\boldsymbol{\sigma}^{*}=\boldsymbol{C}^{*}\left(\boldsymbol{\Gamma}^{*}\right)^{1/2}\Delta
t^{-1/2},$

simulate the processes by using only two stochastic factors and
$\boldsymbol{\sigma}^{*}$ and compute the empirical matrix of log-returns,
obtaining:

$\rho^{*}=\begin{bmatrix}1.00&0.99&0.99&0.97\\\ 0.99&1.00&0.99&0.95\\\
0.99&0.99&1.00&0.95\\\ 0.97&0.95&0.95&1.00\\\ \end{bmatrix},$ (11)

which is extremely close to the one shown in Equation (10).

Figure 10: Possible simulation of the forward prices process $\boldsymbol{F}$
.

If we compute the original covariance matrix $\hat{\boldsymbol{\Sigma}}$ and
$\boldsymbol{\Sigma}^{*}$ given by:

$\boldsymbol{\Sigma}^{*}=\Delta
t\boldsymbol{\sigma}^{*}\left(\boldsymbol{\sigma}^{*}\right)^{T}$

we get:

$\displaystyle\hat{\boldsymbol{\Sigma}}$
$\displaystyle=\begin{bmatrix}1.4859&2.3309&1.6781&1.3756\\\
2.3309&3.6897&2.6580&2.1138\\\ 1.6781&2.6580&1.9197&1.5256\\\
1.3756&2.1138&1.5256&1.3405\\\ \end{bmatrix}\cdot 10^{-4},$
$\displaystyle\boldsymbol{\Sigma}^{*}$
$\displaystyle=\begin{bmatrix}1.4846&2.3302&1.6800&1.3760\\\
2.3302&3.6892&2.6592&2.1140\\\ 1.6800&2.6592&1.9169&1.5250\\\
1.3760&2.1140&1.5250&1.3404\\\ \end{bmatrix}\cdot 10^{-4},$

which confirms that only two Brownian motions are enough to explain the most
part of the variance of log-returns and hence to model the market.

##  4 Calibration

In the previous section, we have shown how to calibrate the model on a toy
example. In this section, we move to a real market application, showing how to
calibrate the model on futures power and gas markets quotations. In most power
markets both peak-load and base-load products are present. The difference
between them is that peak-load contracts delivers electricity only is some
specific hours (typically from 8 to 20) during the working days, whereas a
base-load contract delivers power for all the hours of the delivery periods.
An off-peak contract, which is not typically traded, delivers energy in non
peak-load periods and can be obtained from peak-load and base-load quotations.
Usually traders and risk managers are interested in considering both type of
products. So far, we have not specified the difference between base-load and
peak-load contracts, but we observe we can easily include them. In particular,
in spot simulations both peak-load and base-load products’ volatility might be
included. A possible strategy is to simulate peak-load and base-load spot
prices separately and hence merging together. On the other hand, in many power
futures market peak-load contracts are less liquid that the corresponding
base-load ones and hence “expertise adjustment” has to be done to properly
retrieve peak-load quotations when they are missing. This leads to the
introduction of arbitrary choices and hence to questionable spot simulations
behavior. For these reason and in order to simplify the exposition we focus
only on base-load quotation.

Before calibrating the model we need some data preparation which is crucial to
properly fit the parameters. This part is the most delicate one and great
attention should be payed to it: if the data are not properly prepared,
calibration might leads to wrong parameters estimation and output results
might be misleading. This section is split in two parts. First we show how to
prepare the data, then we use the PCA technique to calibrate the model.

###  4.1 Data preparation

As discussed in Section 3 we have assumed that the volatility
$\boldsymbol{\sigma}(t,T)$ is a function of the time to maturity $T-t$. By
previous discussion, the volatility function must be fitted on fixed maturity
products $F(t,T)$ but, unfortunately, in the market we observe
$F\left(t,T^{s},T^{e}\right)$, namely the price of a futures contract with
maturity $T^{s}$ and delivery period $\left[T^{s},T^{e}\right]$ at tile $t$
(see Table 1).

Trading date | Jan-20 | Feb-20 | Mar-20 | Q2-20 | Q3-20 | Q4-20 | 2021 | 2021
---|---|---|---|---|---|---|---|---
2020-01-02 | 36.05 | 39.76 | 37.15 | 35.50 | 39.05 | 45.30 | 43.85 | 46.55
Table 1: Some of the products $F(t,T^{s},T^{e})$ available on German futures
market for the trading date $4/1/2020$.

Before proceeding further, we need some notation. If today is the 16th
February 2023, we name the product which delivers 1 MWh for all the hours of
the next year $Y1$, namely the closest available calendar traded at the
current date. The same reasoning can be applied to monthly, quarterly and
seasonal products. For example, $Q1$ is the product which delivers in the
following quarter. In our model we are interested in fitting the volatility of
the $M0,M1,\dots,Q1,\dots,Y1,Y2,\dots$. Observe that such products are exactly
the fixed delivery futures contracts $F(t,T)$. Therefore, the first step is to
switch from $F(t,T^{s},T^{e})$ in Table 1 to $F(t,T)$ listed in Table 2.

Observe that if we are on January 2020, $M1$ is February 2020 which
corresponds to the product $F(t,T^{s},T^{e})$ with $T^{s}=1/2/2020$ and
$T^{e}=28/2/2020$ (in red in Tables 1 and 2) . When we move to February 2020
by referring to $M1$ we need to consider the product
$F\left(t,T^{s},T^{e}\right)$ with delivery dates on April 2023 and so on. By
proceeding in this way we construct the time series of the fixed delivery
products $F(t,T)$, namely of $M0,M1,\dots,Q1,\dots,Y1,\dots$. An example of
this structure is shown in Table 2. As stated before, observe that at the end
of March 2020 (green quotation), $Q3$ is the fourth quarter of 2020 whereas at
the beginning of April 2020 the product $Q3$ is the first quarter of the 2021.
This is what is sometimes know as “rolling effect”. The same happens for the
$M1$. In computing the log-returns of fixed delivery products from Table 2
attention should be payed, since when products roll, fake spikes in log-
returns might be created: such spikes as not due to the market structure but
simply on the way we look at it and hence must be removed.

Trading date | $M0$ | $M1$ | $M2$ | $Q1$ | $Q2$ | $Q3$ | $Y1$ | $Y2$ |
---|---|---|---|---|---|---|---|---|---
2020-01-02 | 36.05 | 39.76 | 37.15 | 35.5 | 39.05 | 45.3 | 43.85 | 46.55 |
2020-01-03 | 38.06 | 40.4 | 37.8 | 36.55 | 39.85 | 45.97 | 44.85 | 47.08 |
2020-01-06 | 36.93 | 39.46 | 37.29 | 36 | 39.36 | 45.5 | 44.55 | 46.89 |
… | … | … | … | … | … | … | … | … |
2020-03-30 | 22.49 | 17.6 | 19.08 | 26.74 | 26.84 | 33.74 | 34.72 | 38.3 |
2020-03-31 | 15.74 | 17.06 | 19.79 | 26.74 | 27.63 | 34.45 | 35.65 | 39.05 |
2020-04-01 | 15.74 | 19.04 | 23.24 | 26.74 | 33.76 | 36.24 | 34.95 | 38.56 |
2020-04-02 | 16.39 | 19.09 | 23.26 | 26.86 | 33.98 | 36.57 | 35.35 | 39.04 |
… | … | … | … | … | … | … | … | … |
2023-02-13 | 129.69 | 131.22 | 126.49 | 131.05 | 149.58 | 177.10 | 158.18 | 128.11 |
2023-02-14 | 129.95 | 131.01 | 126.46 | 131.33 | 150.25 | 179.23 | 159.34 | 128.50 |
Table 2: Relative products for the German forward markets.

At this point one might be tempted to compute the log-returns (opportunely
filtered) of the prices in Table 2 and apply the PCA to identify the number of
risky factors. Of course this would be possible if one is interested in the
simulation of the fixed delivery products $M0,M1,\dots,Q1,\dots$ but the
situation is slightly different if one is interested in the simulation of
$F(t,T^{s},T^{e})$, for example Jan-21, on the period starting from the 1st
January 2020 to its maturity which occurs at the 31st December 2020. Observing
Figure 11 and assuming that only $M1$, $M2$, $M3$, $Q1$, $Q2$, $Q3$ and $Y1$
are available, the right volatility of this product is that of the $Y1$ until
the end of March 2020, then we have to consider that of the $Q3$ up to the end
of July and so on until December in which the proper volatility is that of the
$M1$. This observation, which is trivial from a theoretical point of view,
might bring to some thorny problems when we move to the model implementation
since the right volatility must be considered.

$1Jan$$31Mar$$30Jun$$30Sep$$31Oct$$30Nov$$31Dec$$Y+1$$Q+3$$Q+2$$M+3$$M+2$$M+1$
Figure 11: Volatility to use for the simulation of the product Jan-21 over the
period 1/1/2020 to 31/12/2020.

In order to consider the right volatility for each time step of the simulation
many approaches are possible. Observing Figure 11 the simplest thing is to fit
the volatility structure on raw data in Table 2, compute for each time step
how long does it last to the delivery and use the proper corresponding
volatility. For example, if we are simulating Jan-21, in November 2020 we will
use the volatility of the $M2$ product. This approach is the fastest to
implement but it is not free of pitfalls. For example, if we want to introduce
the Brent quotations, we observe that for such a product nor quarter neither
calendar products are quoted, but only months. For this reason, it would not
be very clear how to deal with Brent monthly quotations with with the coarser
granularity ones in power markets insisting on the same delivery period. A
possible solution is to reduce Brent monthly quotations to quarterly or yearly
but in this case we would lose some information.

On the other hand, another approach is to create monthly quotations starting
from the products with a coarser granularity, namely quarter and yearly and
hence fit the volatility parameters on these processed data. In order to do
so, first we have to compute a _flat monthly forward curve_ by obtaining the
value of each monthly futures contract in a market coherent way, guarantying
that no arbitrage opportunities arise. A well known approach for the
construction of a smooth forward curve with a seasonal effect which is
coherent with the marketd observed futures quotations has been proposed by
Benth et al. [9]. A similar procedure can be adopted and simplified by
substituting the smooth curve with a step-wise one. The algorithm, which
reduce to the solution of a linear system, is detailed in the following
session.

####  4.1.1 Construction of the monthly forward curve

In this section we show how to construct the monthly forward curve for a given
market starting from the futures products available on the market. Details can
be found in Benth et al. [9, chapter 7].

Consider the example in Figure 12 where overlapped products are allowed and
suppose we want to define a step-wise forward curve on the intervals
$\left[\tau_{i-1},\tau_{i}\right],\;i=1,\dots,n$ with $n=7$ of the form:

$\epsilon\left(u\right)=\sum_{i=1}^{n}a_{i}\mathbbm{1}_{\left[\tau_{i-1},\tau_{i}\right]}(u),$

where $\left\\{a_{i}\right\\}_{i=1}^{n}$ are the values we have to fit.
Assuming $t_{0}=0$ and denoting the price of the futures product at $t_{0}$ by
$F\left(T_{i}^{s},T_{i}^{e}\right)=F\left(t_{0},T_{i}^{s},T_{i}^{e}\right)$ we
have to guarantee that non arbitrage constraints are satisfied. For each
quoted product $i=1,\dots,n$, we must satisfy the following relation:

$F\left(T_{i}^{s},T_{i}^{e}\right)=\frac{1}{T_{i}^{e}-T_{i}^{s}}\int_{T_{i}^{s}}^{T_{i}^{e}}\epsilon(u)du.$

It is easy to show that finding the values $\left\\{a_{i}\right\\}_{i=1}^{n}$
reduces to the solution of a linear system which can be done numerically in a
very efficient way. Once that the values of $\left\\{a_{i}\right\\}_{i=1}^{n}$
are available $\epsilon(u)$ is determined. For each trading date and for each
market such a monthly forward curve must be computed. If we consider as
trading date $4/1/2020$ the resulting monthly forward curve is shown in Figure
13. Of course, if during the construction of the forward curve some products
are completely overlapped they must be removed preserving those with the finer
granularity.

$F(t,T_{i}^{s},T_{i}^{e})$$T$$T_{1}^{s}$$T_{1}^{e}$$T_{2}^{s}$$T_{2}^{e}$$T_{3}^{s}$$T_{3}^{e}$$T_{4}^{s}$$T_{4}^{e}$$T_{5}^{s}$$T_{5}^{e}$$T_{6}^{s}$$T_{6}^{e}$$T_{7}^{s}$$T_{7}^{e}$$\tau_{0}$$\tau_{1}$$\tau_{2}$$\tau_{3}$$\tau_{4}$$\tau_{5}$$\tau_{6}$$\tau_{7}$
Figure 12: Quoted futures product $F\left(t,T_{i},T_{j}\right)$. Figure 13:
Reconstructed monthly forward curve (black) compared with the quoted forward
price for the trading date 2020-04-01 (green). Observe that non arbitrage
constraints are satisfied on the year 2021.

As useful remark, we observe that no bootstrapping has been done on the
monthly forward curve. The bootstrapping of the forward curve is very common
in energy markets since seasonal effects are present in many countries.
Nevertheless, our final goal is to find the volatility of forward products and
hence a bootstrapping method could introduce some deformation in quotations
which might have an impact on the estimation of the volatility structure
leading to biased results.

Finally, we check if the latter methodology has an impact on the PCA we use to
detect the number of random factors which drive the market. Since we are only
obtaining monthly products form those with a coarser granularity, we are not
expected to introduce more information and variability with respect the case
in which we work with row date of Table 2. In order to verify this claim we
compare the number of factor we need to explain the $95\%$ of the variance
both using raw and monthly data from German power futures market over the
period from $1/1/2020$ to $31/12/2021$. From Figure 14 we observe that in both
cases fourteen principal components are needed to explain around the $95\%$ of
the total variance. In case of monthly data the original total number of
random factors is much larger than that of the raw data, since many quotations
have been replicated on a finer granularity scale. Nevertheless, since many of
them are just replications of exiting ones they are not needed to explain a
significant amount of the whole market variance as confirmed by Figure 15.

Figure 14: Percentage of the total variance explained by principal components
on the raw and on the monthly data. Figure 15: Variance explained from the
principal components of the monthly data. One can easily observe that only the
first principal components are required to explain a sufficiently large amount
of the total variance.

In light of the preceding results we can claim that the calibration of
volatility parameters using the PCA procedure we sketched in Section 4.2, can
be done both on raw data or on processed monthly data. The latter approach is
preferable when commodity with different granularity are present: in this
case, working with products having the same granularity leads to cleaner
analysis and easier coding.

###  4.2 PCA for dimension reduction

In this section we briefly formalize the example in Section 3.4 and we show
how PCA can be used in order to identify a relatively small number of
stochastic factors which drive the whole market. Assume we have $K$ markets
and for each of them we have the same number of contracts with given
maturities and call them $F^{k}(t,T_{m})$, with $k=1,\dots,K$ and
$m=1,\dots,M$. Consider a matrix consisting of $n_{obs}$ rows and
$\tilde{N}=M\cdot K$ columns: each row contains the monthly forward curves
constructed as show in Section 4.1. We compute the matrix of log-returns
$\boldsymbol{X}\in\mathbb{R}^{(n_{obs}-1)\times MK}$. Now we apply the PCA and
we compute the eigenvalues and the associated eigenvectors. We have the
following proposition.

######  Theorem 4.1.

[Johnson and Wichern [29, Result 8.2]] Let
$\boldsymbol{X}=\left(X_{1},\dots,X_{d}\right)$ be a random vector with
covariance matrix $\Sigma$ with eigenvalue-eigenvector pair
$\left(\lambda_{1},v_{1}\right),\dots,\left(\lambda_{d},v_{d}\right)$ where
$\lambda_{1}\geq\lambda_{2}\geq\dots,\geq\lambda_{d}\geq 0$. Let
$Y_{i}=v_{1}\boldsymbol{X}$ for $i=1,\dots,d$ be the principal components.
Then:

$\sigma_{11}+\dots+\sigma_{dd}=\sum_{i=1}^{d}var\left(X_{i}\right)=\lambda_{1}+\dots+\lambda_{d}=\sum_{i=1}^{d}var(Y_{i})$

This result states that the total population variance is given by the sum of
the eigenvalue $\lambda_{i}$, for $i=1,\dots,d$. Hence the percentage of the
variance explained by the principal component $Y_{k}$ is given by:

$\frac{\lambda_{k}}{\sum_{i=1}^{d}\lambda_{i}},\quad k=1,\dots,d.$

The number of principal components $k$ is chosen such that they are enough to
explain a sufficiently large percentage of the variance. Usually this choice
is made heuristically, namely there is not a rigorous way to chose the number
$k$. See Johnson and Wichern [29] for details.

Once that we have chosen how many principal components to consider, say $N$,
we can define the matrix $\boldsymbol{\sigma}^{*}$ as:

$\boldsymbol{\sigma}^{*}=\boldsymbol{C}^{*}\left(\boldsymbol{\Gamma}^{*}\right)^{1/2}\Delta
t^{-1/2},$

where $\boldsymbol{C}^{*}$ is the $\mathbb{R}^{\tilde{N}\times N}$ matrix
consisting on the first $N$ eigenvector associated to the eigenvalues
$\lambda_{1},\dots,\lambda_{N}$, $\boldsymbol{\Gamma}^{*}$ is the diagonal
$\mathbb{R}^{N\times N}$ matrix containing the eigenvalues and $\Delta t$ is
the time interval between quotations (for example $\Delta t=1/252$ in case of
daily quotations). Once we have fitted the matrix $\boldsymbol{\sigma}^{*}$ we
can use Equation (6) to simulate the market by using $N$ independent
stochastic factors.

##  5 Numerical results

In this section we consider a real market application. We calibrate the model
on futures market quotations over the period from $1/1/2020$ to $31/12/2020$.
By considering monthly fixed delivery futures products, obtained as shown in
Section 4.1, in Table 2 for six different markets: four European power futures
markets, Germany (DE), France (F7), Italy (IT) and Switzerland (CH) and two
natural gas markets, PSV and TTF. We consider deliveries up to $M=24$ months
ahead, but larger maturities can be considered if products are available. As
observed in Section 2 we expect a significant level of co-integration between
markets. In Figure 19 we plot log-returns correlation matrices: the level of
historical correlation is high for all maturities and it is higher for the
long-term ones. Moreover, we observe the highest level of correlation across
markets of the same commodity type, namely power and natural gas. Furthermore,
the correlation between power and natural gas is significant too, since the
natural gas is commonly used as fuel for electricity production in many
European countries.

Starting from this quotations we compute the log-returns matrix
$\boldsymbol{X}$. As discussed in the previous section, “fake” spikes in log-
returns might be created and they must be removed. Moreover, since we are
assuming that log-returns are Gaussian, we filter out the outliers by removing
everything which lies more that three standard deviation away from the mean.
This is a very rough way of filtering outliers out but it is a method commonly
used by practitioners. Otherwise a similar approach to the one proposed by
Cartea and Figueroa [19] can be adopted. Once matrix $\boldsymbol{X}$ has been
prepared, we perform a PCA analysis and we select a sufficient bulk of
stochastic factors which are enough to explain a sufficiently large part of
the total market variance. As observed by Feron and Gruet [23], Geman [25] and
Koekebakker and Fridthjof [31] 5-10 factors in power markets are enough to
explain a large part of the variance, whereas for other markets, such as
metals, 4 are enough. In this example we start from 144 stochastic factors
and, after a PCA, analysis we get that 10 factors are enough to explain the
90% of the total market variance.

Figure 16: Simulated, market and theoretical annualized volatility.

Following the procedure we illustrated in Section 3.4, by PCA we calibrate the
model by fitting the matrix $\boldsymbol{\sigma}^{*}$. Once this has been
done, we are ready to simulate the structure of the forward curve. Considering
as $t_{0}=4/1/2021$, by using Equation (6) we simulate futures price
$F^{k}(t,T^{s},T^{e})$ with delivery period $T^{s}=1/1/2022$,
$T^{e}=31/12/2022$ form $t_{0}$ up to the $31^{st}$ December 2021. The result
is shown in Figure 17. The dynamic of the prices exhibits an extremely high
level of dependence: in particular, both power and natural gas prices tends to
rise or fall together over the whole simulation period.

Figure 17: Sample paths for the product calendar 2022 with delivery date from
2022-01-01 to 2022-12-31.

In Figure 18 we show some possible path realizations on a time interval
$\left[t_{0},T\right]$ , together with the fifth and the ninety-fifty
percentile, of the single product German power futures calendar
$F(t,T^{s},T^{e})$ with delivery dates $T^{s}=1/1/2022$ and
$T^{e}=31/12/2022$. Since the dynamic introduced by Equation (4) is diffusive
we get, as expected, a time increasing variance of the process at time, for
$t\in\left[t_{0},T\right]$.

Figure 18: Sample paths of the product German power futures calendar 2022. In
red the fifth and the ninety-fifth percentile.

Furthermore, in Figure 20 we show the log-returns matrices correlation
computed on simulated fixed delivery products
$F^{k}(t,T_{m}),\;k=1,\dots,K\;m\in\left[1,M\right]$ on a given time interval.
A comparison between Figure 19 and Figure 20 shows that the correlation
structure is properly replicated from the proposed model. Of course, since we
used only 10 stochastic instead of the 144 original ones, the correlation
surface computed from the simulation appears to be less varied than the
original one. In particular, the correlation surface computed from the
simulations is flat for long maturities. This happens because PCA selects a
single stochastic factor to move all the long term structure. This is in
accordance with empirical evidence since products with long time to maturity
are strongly correlated and tend to oscillate in the same way.

Figure 19: Historical daily log-returns correlation surfaces for different
commodities.

Figure 20: Simulated daily log-returns correlation surfaces for different
commodities.

In many cases, practitioners are interested in computing the value at risk of
a given portfolio. In order to achieve this task, many approaches, historical,
parametric and Monte Carlo among the others, are available. If a Monte Carlo
approach is chosen the simulation of the forward curve must be performed. In
order to do so, we use Equation (7) to simulate the forward curve after a few
days, typically one or two. Results are shown in Figure 21. In particular, we
observe that the expected value of the simulated forward curves $\mu$,
converges at today’s forward curve $F(t_{0},T)$ as one should expect by
construction. The fifth and the ninety-fifth percentiles in red, give an idea
of the amplitude of the simulations: this is a useful visual check in order to
show if the simulation routine performed well and it is largely used by
practitioners.

Figure 21: Simulated two days forward curves: $n_{sim}=400$.

As final issue, we investigate the spot prices produced by the HJM framework.
In this step, for the sake of conciseness, we focus only of two markets:
German and Italian power spot markets. In Figure 22 we show a single
realization of the stochastic process
$S^{k}=\left\\{S^{k}(t);t\in\left[t_{0},T\right]\right\\}$ for the two
different commodities, according to the dynamic presented by Equation (9). As
we can observed the spot prices tends to move together since they are driven
by the same stochastic factors. In this case the correlation level in daily
spot log-returns is high, approximately $\rho=0.98$ and this is confirmed by
Figure 23, where we displayed the contour plot of a bi-variate Gaussian
distribution fitted on daily log-returns.

It is worth remembering that many European power spot markets are coupled
together: market coupling optimizes the allocation of cross-border capacities
between countries. One of the possible effects on the prices dynamic is that
at the same hour the electricity price in two different countries is the same.
The proposed model does not take in account this behavior: in order to include
such effect a fundamental component must introduced, as proposed by Carmona
and Coulon [18] and Kiesel and Kusterman [30]. Another limitation of the
approach we proposed is that the correlation we introduce in the spot
simulations is the one we derive from futures prices which is typically higher
than the one we get historically in the spot market. Despite this limitation,
the methodology the methodology we introduced is a good way to produce a
coherent framework both for spot and futures energy markets.

Finally we stress out that, since Monte Carlo simulations are available both
for power and futures products derivatives, pricing of complex energy
financial claims, such as virtual power plants or storage might be easily
performed following the algorithms presented by Tseng and Barz [50] and
Boogert and de Jong [12] respectively.

(a) Spot prices simulations.

(b) Zoom spot prices simulations.

Figure 22: German and Italian power spot prices simulations. Figure 23:
Contour plot of a bi-variate normal distribution fitted on daily log-returns
of German and Italian power spot prices simulation. Correlation in log-returns
$\rho=0.98$. Figure 24: Sample paths for the German power spot price.

##  6 Conclusions

In this article we discussed in detail the implementation and a possible
application of the Heat-Jarrow-Morton framework to energy markets from a very
pragmatic point of view. In particular we have focused on the European power
and natural gas markets. We introduced a Black-style dynamic for the so called
_fixed-maturity_ products and we derived the ones of spot and futures prices.
By showing that all the power and natural gas futures markets we considered
are strictly dependent, we selected by the PCA algorithm only few stochastic
factors which explain a large part of the variance. Furthermore, following
Böerger et al. [17] we have shown how a closed form solution for European
vanilla option is available. Finally, we have applied the model to real market
data from European power and natural gas markets. We discussed the daily log-
returns correlation structures and we have shown that the models fits the
market narrowly. Moreover, we analyzed the futures and spot simulations in
output from the model.

HJM model is easy to implement, calibration step is not difficult, due to the
hypothesis of normality of log-returns and simulations are fast to perform
without introducing a discretization error, since an exact solution for the
stochastic differential equation is available.

Unfortunately, the model presents some drawbacks which must be considered when
the output are used for risk-metric computation or for derivative pricing.
First of all, log-returns are normally distributed and this is not the case in
almost all financial markets: jumps, volatility smiles and clustering are
often present and are not considered by the proposed framework. In order to
consider such stylized facts one can introduce stochastic volatility, jumps in
price dynamics maybe using subordination techniques. On the other hand, even
if such models work great in uni-variate setting, they are too hard to
calibrate in a multivariate framework, especially when the number of
underlying assets is higher than three. In this latter case, a Gaussian
approach could be preferred in term of computational al calibration
complexity. On the other hand, if one needs to focus on the pricing of a
derivative written o a single underlying asset, more complex models based on
Lévy processes or on stochastic volatility can be used in order to include
many market stylized facts.

Other drawbacks follow from the fact that the calibration has been performed
on historical futures prices. First of all, it is not guaranteed that the
model replicates the options quoted in the market: in order to guarantee that,
a risk-neutral calibration should be performed. Unfortunately, option power
markets are not very liquid and hence, in many situations, an historical
calibration is the only possible way. The second arising problem is that,
since we did not considered spot quotations, the correlation and the
volatility in simulated log-returns directly inherit from the forward one
instead of from the spot one.

On the other hand, the proposed framework appears to be one of the easiest
technique to obtain dependent futures and spot prices simulation in a multi-
commodity setting in a market coherent way and, for this reason, it is widely
used by practitioners. Nevertheless, it would be worth to include jumps in
prices’ dynamic or to consider a stochastic volatility approach in the multi-
dimensional framework in order to better model the prices’ dynamic, by
preserving both mathematical and numerical tractability. This could be the
direction of possible future researches.

## References

  * Applebaum [2009] D. Applebaum.  _Lévy Processes and Stochastic Calculus_.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2 edition, 2009.  doi: 10.1017/CBO9780511809781.
  * Asche et al. [2006] F. Asche, P. Osmundsen, and M. Sandsmark.  The UK Market for Natural Gas, Oil and Electricity: Are the Prices Decoupled?  _The Energy Journal_ , 27(2):27–40, 2006.  ISSN 01956574, 19449089.  URL http://www.jstor.org/stable/23297017.
  * Bachmeier and Griffin [2006] Lance J. Bachmeier and James M. Griffin.  Testing for market integration crude oil, coal, and natural gas.  _The Energy Journal_ , 27(2):55–71, 2006.  ISSN 01956574, 19449089.  URL http://www.jstor.org/stable/23297019.
  * Ballotta and Bonfiglioli [2013] L. Ballotta and E. Bonfiglioli.  Multivariate Asset Models Using Lévy Processes and Applications.  _The European Journal of Finance_ , 13(22):1320–1350, 2013.
  * Barndorff-Nielsen [1998] O.E. Barndorff-Nielsen.  Processes of Normal Inverse Gaussian Type.  _Finance and Stochastics_ , 2(1):41–68, 1998\.
  * Bates [1996] D. S. Bates.  Jumps and Stochastic Volatility: Exchange Rate Processes Implicit in Deutsche Mark Options.  _Review of Financial Studies_ , 9(1):69–107, 1996\.
  * Bencivenga et al. [2010] C. Bencivenga, G. Sargenti, and R. L. D’Ecclesia.  Energy Markets: Crucial Relationship Between Prices.  In Marco Corazza and Claudio Pizzi, editors, _Mathematical and Statistical Methods for Actuarial Sciences and Finance_ , pages 23–32, Milano, 2010. Springer Milan.  ISBN 978-88-470-1481-7.
  * Benth and Saltyte-Benth [2007] F.E. Benth and J. Saltyte-Benth.  The Volatility of Temperature and Pricing of Weather Derivatives.  _Quantitative Finance_ , 7(5):553–561, 2007.  doi: 10.1080/14697680601155334.
  * Benth et al. [2008] F.E. Benth, J.S. Benth, and S. Koekebakker.  _Stochastic Modeling of Electricity and Related Markets_.  Number 6811 in World Scientific Books. World Scientific Publishing Co. Pte. Ltd., July 2008.  ISBN ARRAY(0x50724998).  URL https://ideas.repec.org/b/wsi/wsbook/6811.html.
  * Benth et al. [2017] F.E. Benth, M. Piccirilli, and T. Vargiolu.  Additive Energy Forward Curves in a Heath-Jarrow-Morton Framework.  Papers 1709.03310, arXiv.org, September 2017.  URL https://ideas.repec.org/p/arx/papers/1709.03310.html.
  * Black and Scholes [1973] F. Black and M. Scholes.  The Pricing of Options and Corporate Liabilities.  _Journal of Political Economy_ , 81(3):637–654, 1973.
  * Boogert and de Jong [2008] A. Boogert and C. de Jong.  Gas Storage Valuation Using a Monte Carlo Method.  _Journal of Derivatives_ , 15:81–91, 2008.
  * Bosco et al. [2010] B. Bosco, L. Parisio, M. Pellegatti, and F. Baldi.  Long-run Relations in European Electricity Prices.  _Journal of Applied Econometrics_ , 25(5):805–832, 2010.  ISSN 08837252, 10991255.  URL http://www.jstor.org/stable/40865175.
  * Brigo et al. [2007] D. Brigo, A. Dalessandro, Neugebauer M., and F. Triki.  A Stochastic Process Toolkit for Risk Management.  Technical report, 2007.
  * Broszkiewicz-Suwaj and Aleksander [2005] E. Broszkiewicz-Suwaj and W. Aleksander.  Calibration of the Multifactor HJM Model for Energy Market.  HSC Research Reports HSC/05/03, Hugo Steinhaus Center, Wroclaw University of Technology, 2005.  URL https://EconPapers.repec.org/RePEc:wuu:wpaper:hsc0503.
  * Buchmann et al. [2019] B. Buchmann, K. Lu, and D. Madan.  Calibration for Weak Variance-Alpha-Gamma Processes.  _Methodology and Computing in Applied Probability_ , 21(4), 2019.  doi: 10.1007/s11009-018-9655-y.
  * Böerger et al. [2009] Reik H. Böerger, Rüdiger Kiesel, and Gero Schindlmayr.  A Two-Factor Model for the Electricity Forward Market.  _Quantitative Finance_ , 9(3):279–287, 2009.  doi: https://doi.org/10.1080/14697680802126530.
  * Carmona and Coulon [2014] René Carmona and Michael Coulon.  _A Survey of Commodity Markets and Structural Models for Electricity Prices_ , pages 41–83.  Springer New York, New York, NY, 2014.  ISBN 978-1-4614-7248-3.  doi: 10.1007/978-1-4614-7248-3˙2.  URL https://doi.org/10.1007/978-1-4614-7248-3_2.
  * Cartea and Figueroa [2005] A. Cartea and M. Figueroa.  Pricing in Electricity Markets: a Mean Reverting Jump Diffusion Model with Seasonality.  _Applied Mathematical Finance, No. 4, December 2005_ , 12(4):313–335, 2005.
  * Cont and Tankov [2003] R. Cont and P. Tankov.  _Financial Modelling with Jump Processes_.  Chapman and Hall, 2003.
  * Edoli et al. [2013] E. Edoli, D. Tasinato, and T. Vargiolu.  Calibration of a multifactor model for the forward markets of several commodities.  _Optimization_ , 62(11):1553–1574, 2013.  doi: 10.1080/02331934.2013.854786.  URL https://doi.org/10.1080/02331934.2013.854786.
  * Emery and Liu [2002] G. W. Emery and Q. Liu.  An Analysis of the Relationship Between Electricity and Natural-Gas Futures Prices.  _Journal of Futures Markets_ , 22(2):95–122, 2002\.  doi: https://doi.org/10.1002/fut.2209.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fut.2209.
  * Feron and Gruet [2020] O. Feron and P. Gruet.  Estimation of the Number of Factors in a Multi-Factorial Heath-Jarrow-Morton model in Electricity Markets.  Working Papers hal-02880824, HAL, June 2020.  URL https://ideas.repec.org/p/hal/wpaper/hal-02880824.html.
  * Gardini et al. [2021] M. Gardini, P. Sabino, and E. Sasso.  Correlating Lévy Processes with Self-Decomposability: Applications to Energy Markets.  _Decision in Economics and Finance_ , 2021.
  * Geman [2005] H. Geman.  _Commodities and Commodity Derivatives. Modeling and Pricing for Agriculturals, Metals and Energy_.  John Wiley & Sons, 2005.  URL https://shs.hal.science/halshs-00144182.
  * Heath et al. [1992] D. Heath, R. Jarrow, and A. Morton.  Bond Pricing and the Term Structure of Interest Rates: A New Methodology for Contingent Claims Valuation.  _Econometrica_ , 60(1):77–105, 1992.
  * Heston [1993] S. L. Heston.  A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options.  _The Review of Financial Studies_ , 6(2):327–343, 1993.
  * Hinderks et al. [2020] W.J. Hinderks, R. Korn, and A. Wagner.  A Structural Heath–Jarrow–Morton Framework for Consistent Intraday Spot and Futures Electricity Prices.  _Quantitative Finance_ , 20(3):347–357, 2020\.  doi: 10.1080/14697688.2019.1687927.  URL https://doi.org/10.1080/14697688.2019.1687927.
  * Johnson and Wichern [2002] R.A. Johnson and D.W Wichern.  _Applied multivariate statistical analysis_.  Prentice Hall, Upper Saddle River, NJ, 5. ed edition, 2002.  ISBN 0130925535.  URL http://gso.gbv.de/DB=2.1/CMD?ACT=SRCHA&SRT=YOP&IKT=1016&TRM=ppn+330798693&sourceid=fbw_bibsonomy.
  * Kiesel and Kusterman [2016] R. Kiesel and M. Kusterman.  Structural Models for Coupled Electricity Markets.  _Journal of Commodity Markets_ , 3(1):16–38, 2016\.  ISSN 2405-8513.  doi: https://doi.org/10.1016/j.jcomm.2016.07.007.  URL https://www.sciencedirect.com/science/article/pii/S2405851315300933.
  * Koekebakker and Fridthjof [2005] S. Koekebakker and O. Fridthjof.  Forward Curve Dynamics in the Nordic Electricity Market.  _Managerial Finance_ , 31:73–94, 2005.
  * Lucia and Schwartz [2002] J. Lucia and E. Schwartz.  Electricity Prices and Power Derivatives: Evidence from the Nordic Power Exchange.  _Review of Derivatives Research_ , 5(1):5–50, Jan 2002.  ISSN 1573-7144.  doi: 10.1023/A:1013846631785.  URL https://doi.org/10.1023/A:1013846631785.
  * Luciano and Semeraro [2010] E. Luciano and P. Semeraro.  Multivariate Time Changes for Lévy Asset Models: Characterization and Calibration.  _Journal of Computational and Applied Mathematics_ , 233(1):1937–1953, 2010.
  * Madan et al. [1999] D. Madan, P. Carr, M. Stanley, and E. Chang.  The Variance Gamma Process and Option Pricing.  _Review of Finance_ , 2, 11 1999.  doi: 10.1023/A:1009703431535.
  * Merton [1976] R.C. Merton.  Options Pricing when Underlying Shocks are Discontinuous.  _Journal of Financial Economics_ , 3:125–144, 1976.
  * Mjelde and Bessler [2009] J. W. Mjelde and D.A. Bessler.  Market Integration Among Electricity Markets and Their Major Fuel Source Markets.  _Energy Economics_ , 31(3):482–491, 2009.  ISSN 0140-9883.  doi: https://doi.org/10.1016/j.eneco.2009.02.002.  URL https://www.sciencedirect.com/science/article/pii/S0140988309000310.
  * Panagiotidis and Rutledge [2007] T. Panagiotidis and E. Rutledge.  Oil and Gas Markets in the UK: Evidence from a Cointegrating Approach.  _Energy economics_ , 29(2):329–347, 2007.  ISSN 0140-9883.
  * Panov and Samarin [2019] V. Panov and E. Samarin.  Multivariate Asset-Pricing Model Based on Subordinated Stable Processes.  _Applied Stochastic Models in Business and Industry_ , 35(4):1060–1076, 2019.
  * Park et al. [2006] H. Park, J. W. Mjelde, and D. A. Bessler.  Price Dynamics Among U.S. Electricity Spot Markets.  _Energy Economics_ , 28(1):81–101, 2006.  ISSN 0140-9883.  doi: https://doi.org/10.1016/j.eneco.2005.09.009.  URL https://www.sciencedirect.com/science/article/pii/S0140988305000903.
  * [40] N. Cufaro Petroni and P. Sabino.  Cointegrating Jumps: an Application to Energy Facilities.  accepted presentation at the $10^{t}h$ Energy and Finance Conference, London, 9-11 September 2015.
  * Piccirilli et al. [2021] M. Piccirilli, M. D. Schmeck, and T. Vargiolu.  Capturing the Power Options Smile by an Additive Two-Factor Model for Overlapping Futures Prices.  _Energy Economics_ , 95:105006, 2021.  ISSN 0140-9883.  doi: https://doi.org/10.1016/j.eneco.2020.105006.  URL https://www.sciencedirect.com/science/article/pii/S0140988320303467.
  * Saifert and Uhrig-Homburg [2006] J. Saifert and M. Uhrig-Homburg.  Modelling Jumps in Electricity Prices: Theory and Empirical Evidence.  2006\.
  * [43] P. A. Samuelson.  _Proof that Properly Anticipated Prices Fluctuate Randomly_ , chapter Chapter 2, pages 25–38.  doi: 10.1142/9789814566926˙0002.  URL https://www.worldscientific.com/doi/abs/10.1142/9789814566926_0002.
  * Sato [1999] K. Sato.  _Lévy Processes and Infinitely Divisible Distributions_.  Cambridge U.P., Cambridge, 1999.
  * Schoutens [2003] W. Schoutens.  _Lévy Processes in Finance: Pricing Financial Derivatives_.  John Wiley and Sons Inc, 2003.
  * Schwartz [1997] E. Schwartz.  The Stochastic Behavior of Commodity Prices: Implications for Valuation and Hedging.  _The Journal of Finance_ , 52(3):923–973, 1997\.  doi: https://doi.org/10.1111/j.1540-6261.1997.tb02721.x.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1540-6261.1997.tb02721.x.
  * Schwartz and Smith [2000] P. Schwartz and J.E. Smith.  Short-term Variations and Long-term Dynamics in Commodity Prices.  _Management Science_ , 46(7):893–911, 2000.
  * Sclavounos and Ellefsen [2009] P.D. Sclavounos and P. E. Ellefsen.  Multi-Factor Model of Correlated Commodity Forward Curves for Crude Oil and Shipping Markets.  2009\.
  * Seydel [2004] R. Seydel.  _Tools for Computational Finance_.  Universitext (1979). Springer, 2004.  ISBN 9783540406044.
  * Tseng and Barz [2000] C. Tseng and G. Barz.  Short-Term Generation Asset Valuation: A Real Options Approach.  _Operations Research_ , 50(2):297–310, 2000.