# From the Samuelson Volatility Effect to a Samuelson Correlation Effect: Evidence from Crude Oil Calendar Spread Options

## Abstract

We introduce a multi-factor stochastic volatility model based on the CIR/Heston stochastic volatility process. In order to capture the Samuelson effect displayed by commodity futures contracts, we add expiry-dependent exponential damping factors to their volatility coefficients. The pricing of single underlying European options on futures contracts is straightforward and can incorporate the volatility smile or skew observed in the market. We calculate the joint characteristic function of two futures contracts in the model in analytic form and use the one-dimensional Fourier inversion method of Caldana and Fusai (2013) to price calendar spread options. The model leads to stochastic correlation between the returns of two futures contracts. We illustrate the distribution of this correlation in an example. We then propose analytical expressions to obtain the copula and copula density directly from the joint characteristic function of a pair of futures. These expressions are convenient to analyze the term-structure of dependence between the two futures produced by the model. In an empirical application we calibrate the proposed model to volatility surfaces of vanilla options on WTI. In this application we provide evidence that the model is able to produce the desired stylized facts in terms of volatility and dependence. In a separate appendix, we give guidance for the implementation of the proposed model and the Fourier inversion results by means of one and two-dimensional FFT methods.

Keywords: Commodities Crude Oil Futures Curve Stochastic Volatility Multi-Factor Model Characteristic Function Fourier Transform Calendar Spread Option

JEL: C63 C52 G13

## 1 Introduction

Crude oil is by far the world’s most actively traded commodity. It is usually traded on exchanges in the form of futures contracts. The two most important benchmark crudes are West Texas Intermediate (WTI), traded on the NYMEX, and Brent, traded on the ICE. Recently, the Dubai Mercantile Exchange’s (DME) Oman contract has been attracting investors looking for a Middle Eastern sour crude oil benchmark. In the S&P Goldman Sachs Commodity Index, WTI has a weight of and Brent a weight of , for a combined total of almost half the index. Another widely quoted index, Jim Rogers’ RICI, has weights of for WTI and for Brent. The crude oil derivatives market is also the most liquid commodity derivatives market. Popular products are European, American, Asian, and calendar spread options on futures contracts.

An important empirical feature of crude oil markets is the absence of seasonality, which is in marked contrast to, say, agricultural commodities markets. A second empirical feature is stochastic volatility of futures contracts, which is clearly reflected in the oil volatility index (OVX), or “Oil VIX”, introduced on the CBOE in July 2008. A third feature is known as the Samuelson effect (Samuelson, 1965; Bessembinder et al., 1996), i.e. the empirical observation that a given futures contract increases in volatility as it approaches its maturity date. Finally, European and American options on futures (usually specified to expire just a couple of days before the underlying futures contract itself) tend to show a more or less strongly pronounced volatility smile, the shape of which depends on the option’s maturity.

European and American options depend on the evolution of just one underlying futures contract. In contrast to these, calendar spread options have a payoff that is calculated from the difference of two futures contracts with different maturities. Therefore, a mathematical analysis and evaluation of calendar spread options must be carried out in a framework that models the joint stochastic behaviour of several futures contracts.

In this article, we propose a multi-factor stochastic volatility model for the crude oil futures curve. Like the popular Clewlow and Strickland (1999a, b) models, the model is futures-based, not spot-based, which means it can trivially match any given futures curve by accordingly specifying the futures’ initial values without “using up” any of the other model parameters. The variance processes are based on the Cox et al. (1985) and Heston (1993) stochastic variance process. However, in order to capture the Samuelson effect, we add expiry dependent exponential damping factors. As in the Heston (1993) model, futures returns and variances are correlated, so that volatility smiles of American and European options observed in the market can be closely matched. The instantaneous correlation of the returns of two futures contracts is also stochastic in our multi-factor model, since it is calculated from the stochastic variances.

Our first result is the calculation of the joint characteristic function of the log-returns of two futures contracts in analytic form. Using this function, calendar spread option prices can be obtained via -dimensional Fourier integration as shown by Caldana and Fusai (2013) or the -dimensional Fast Fourier Transform (FFT) algorithm of Hurd and Zhou (2010). The fast speed of these algorithms is of great importance when calibrating the model to these products.

Our second result is to give analytical formulas to recover the dependence structure of two futures prices from the joint characteristic function of the model. The proposed expressions give the dependence structure as the copula function of the two prices or as its copula density. In many studies, the measure chosen to describe dependence is Pearson’s rho, which, however, also depends on the marginal distributions. In order to completely insulate our analysis from the influence of the marginal distributions, we carry it out via the copula function produced by the joint characteristic function. Once we have the copula and its copula density it becomes possible to compute various dependence and concordance measures, such as Spearman’s rho and Kendall’s tau, for two futures contracts.

Using these mathematical tools, we carry out an analysis of the dependence of the returns of two futures contracts. We observe that, for a fixed time-horizon, these returns become less dependent as the maturity of the second underlying futures contract increases and moves away from that of the first underlying contract. In analogy to the classic Samuelson volatility effect, we call this effect the Samuelson correlation effect.

Copula functions can also be used to give a rigorous definition of the implied correlation of calendar spread options. The traditional definition assumes a bivariate Black-Scholes-Merton model for the two underlyings, which assumes in particular that the marginal distributions are log-normal. In contrast, here, following Tavin (2014), and using the actual marginal distributions of the model, for a given calendar spread option price we define the implied correlation as the value of the correlation parameter in the bivariate Gaussian copula that reproduces this price. Note that implied correlation depends both on the strike and the maturity of the option, phenomena usually referred to as correlation smile/skew/frown and correlation term structure.

In an empirical section we calibrate the two-factor version of our model to market data from three different dates. We show that the model can fit European and American option prices as well as calendar spread option prices very closely. The model could therefore be used by a price maker in a crude oil market to provide consistent and arbitrage free prices to other market participants.

We conclude this introduction with a survey of previous literature. A detailed exposition of commodity models is given by Clark (2014). One of the most important and still widely used models is the Black (1976) futures model, which is set in the Black-Scholes-Merton framework. Contracts with different maturities can have different volatilities in this model, but for each contract the volatility is constant. Therefore, Black’s model doesn’t capture the Samuelson effect. European option prices in this model are given by the Black-Scholes-Merton formula, and consequently there is no volatility smile for options with different strikes. Finally, all contracts are perfectly correlated in this model, since they are driven by the same Brownian motion.

Clewlow and Strickland (1999a, b) propose one-factor and multi-factor models of the entire futures curve with deterministic time-dependent volatility functions. A popular specification for these functions is with exponential damping factors. Since this specification still leads to log-normally distributed futures prices, there is no volatility smile or skew in this model. In the one-factor model, the instantaneous returns of contracts with different maturities are perfectly correlated; in the multi-factor model, however, these returns are not perfectly, but deterministically correlated.

Stochastic volatility models have been proposed by Scott (1987, 1997), Wiggins (1987), Hull and White (1987), Stein and Stein (1991), Heston (1993), Bakshi et al. (1997), and Schoebel and Zhu (1999), among others. Extending the Heston (1993) model to multiple factors, Christoffersen et al. (2009) show that under certain independence assumptions it is straightforward to obtain the characteristic function in closed form and calculate European option prices using the Fourier transform. An important aspect they then proceed to study is the stochastic correlation between the stock return and variance implied by the model. Duffie et al. (2000) have studied a very general class of jump-diffusions. The model presented in this paper fits into this framework (in the version extended to time-dependent parameters).

Trolle and Schwartz (2009) introduce a very general two-factor spot based model, with, in addition, two stochastic volatility factors as well as two stochastic factors for the forward cost-of-carry. The variance processes are extensions of the CIR/Heston process to a more general mean-reversion specification. They also give the dynamics of their model in terms of the futures curve. The main focus of their study is on unspanned stochastic volatility of single-underlying options on futures contracts.

Spread options have been well studied in a two-factor Black-Scholes-Merton framework. Margrabe (1978) gives an exact formula when the strike equals zero, and Kirk (1995), Carmona and Durrleman (2003), Bjerksund and Stensland (2011) and Venkatramanan and Alexander (2011) give approximation formulas for any .

Caldana and Fusai (2013) have recently proposed a very fast one-dimensional Fourier method that extends the approximation given by Bjerksund and Stensland (2011) for the Black-Scholes-Merton model to any model for which the joint characteristic function is known. As in Bjerksund and Stensland (2011), the method provides a lower bound for the spread option price, but in practice the bound seems to be so close to the actual option price that it can be used as the price itself.

Carr and Madan (1999) show how the Fast Fourier Transform (FFT) can be used to price European options with different strikes in one step. Dempster and Hong (2002) and Hurd and Zhou (2010) apply the two-dimensional FFT to the pricing of spread options. Hurd and Zhou’s method returns spread option prices at many different strikes (after a re-scaling and interpolation step), in analogy to Carr and Madan (1999), in one inversion step.

The rest of the paper proceeds as follows. In Section 2 we define the proposed model and provide the associated joint characteristic function. Section 3 deals with spread options and the structure of dependence produced by the model. Section 4 presents an empirical analysis based on different market situations. Section 5 concludes.

## 2 A Model with Stochastic Volatility for Crude Oil Futures

### 2.1 The Financial Framework and the Model

We begin by giving a mathematical description of our model under the risk-neutral measure . Let be an integer, and let be Brownian motions under . Let be the maturity of a given futures contract. The futures price at time , is assumed to follow the stochastic differential equation (SDE)

(1) |

The processes are CIR/Heston square-root stochastic variance processes assumed to follow the SDE

(2) |

For the correlations, we assume

(3) |

and that otherwise the Brownian motions are independent of each other. As we will see, this assumption has as a consequence that the characteristic function factors into separate expectations. Note also that the identity together with Sylvester’s criterion can be used to show that the correlation matrix determined by (3) is indeed positive definite for any choice of the parameters .

For fixed , the futures log-price follows the SDE

(4) |

Integrating (4) from time up to a time , gives

(5) |

We define the log-return between times and of a futures contract with maturity as

In the following, the joint characteristic function of two log-returns will play an important role. For , is given by

(6) |

The joint characteristic function of the futures log-prices is then given by

(7) |

Note that futures prices in our model are not mean-reverting, and that the log-price at time and the log-return are independent random variables.

In the following proposition, we show how the joint characteristic function , and therefore also the single characteristic function , is given by a system of two ordinary differential equations (ODE).

###### Proposition 2.1

The joint characteristic function at time for the log-returns of two futures contracts with maturities is given by

where

and the functions and satisfy the two differential equations

with

The single characteristic function at time for the log-return of a futures contract with maturity is given by setting in the joint characteristic function.

The statement regarding the single characteristic function immediately follows from the definition of the joint characteristic function. The joint characteristic function is calculated in appendix A.

In the next proposition, we show how this ODE system can be solved analytically. A closed form expression for is found thanks to a computer algebra software and is then proportional to the integral of on .

###### Proposition 2.2

Dropping the references to , the function is given in closed form as

with and constants with respect to , defined as

The functions and are the confluent hypergeometric functions.

and are usually referred to as Kummer’s functions as they solve Kummer’s equation (Kummer, 1836; Tricomi, 1955). The function is also known as , and the function as Tricomi’s function. Given, , Kummer’s equation is

(8) |

A way to obtain is by means of a series expansion

(9) |

And is obtained from as

(10) |

where denotes the Gamma function extended to the complex plane. These results and additional properties of Kummer’s functions (e.g. integral representations) can be found in Chap. 13 of Abramovitz and Stegun (1972). A detailed analysis of how to implement Kummer’s functions is given by Pearson (2009). A suitable way to implement the complex Gamma function is the Lanczos (1964) approximation.

As has already been mentioned, the model introduced above is an extension of the one-factor Clewlow and Strickland (1999a) and multi-factor Clewlow and Strickland (1999b) models to stochastic volatility. Since these models are useful benchmarks, we give a description of them and calculate their joint characteristic function. In the risk-neutral measure , the futures price is modelled with deterministic time-dependent volatility functions :

(11) |

where are independent Brownian motions. A popular specification for the volatility functions is

(12) |

for fixed parameters , so that the volatility of a contract a long time away from its maturity is damped by the exponential factor(s).

The marginal and joint distributions of futures prices are log-normal in the Clewlow-Strickland models. Nevertheless, it can be very useful to know the single and joint characteristic functions as well for testing and benchmarking purposes. The next proposition gives closed-form solutions for them.

###### Proposition 2.3

In the Clewlow and Strickland model defined by (11) and (12), the joint characteristic function at time for the log-returns of two futures contracts with maturities is given by

The single characteristic function at time for the log-return of a futures contract with maturity is given by setting in the joint characteristic function.

We prove this result in appendix A.

This result can also be used to add non-stochastic volatility factors to the model by multiplying the joint characteristic function of Proposition 2.1 with one or more factors from Proposition 2.3. Since each “Clewlow-Strickland” factor depends on only two parameters and , it does not add a significant burden to the calibration to market data, while allowing for increased flexibility when fitting the model to the observed volatility term structure.

### 2.2 Pricing Vanilla Options

European options on futures contracts can be priced using the Fourier inversion technique as described in Heston (1993) and Bakshi and Madan (2000), or the FFT algorithm of Carr and Madan (1999). Alternatively, they can be priced by Monte Carlo simulation using discretizations of (1) (Euler scheme) or (4) (Log-Euler scheme) and of (2).

Let denote the strike and the maturity of a European call option on a futures contract with maturity , and let the single characteristic function of the futures log-price be given by . In the general formulation of Bakshi and Madan (2000), the numbers

(13) | ||||

(14) |

represent the probabilities of finishing in-the-money at time in case the futures itself or a risk-free bond is used as numéraire, respectively. The price of a European call option is then obtained with the formula

(15) |

European put options can be priced via put-call parity .

American call and put options can be evaluated via Monte-Carlo simulation using the method of Longstaff and Schwartz (2001). Alternatively, the early exercise premium can be approximated with the formula of Barone-Adesi and Whaley (1987). Trolle and Schwartz (2009) (Appendix B) address the issue of estimating European prices from American prices.

A typical WTI volatility surface displays high implied volatilities at the short end and low implied volatilities at the long end. This is in line with the Samuelson effect. Furthermore, there is usually a strongly pronounced smile at the short end, and a weak smile at the long end.

### 2.3 Stochastic Correlation in the Multi-Factor Model

We will show in this section that if we specify our model with two or more volatility factors, then the returns of two given futures contracts are stochastically correlated, which is a realistic and important feature.

Define

(16) |

Then the instantaneous correlation at time is given by:

(17) |

Let us begin with an examination of the -factor model, in which futures returns follow the SDE

(18) |

and the variance process follows the SDE

(19) |

The correlation is given by . Inserting (18) into (16) gives for the instantaneous covariance

(20) |

Cox et al. (1985) show that the random variable follows a non-central -distribution. It is easy to see that the instantaneous correlation (17) is always equal to one in the -factor model:

(21) |

Finally, the terminal covariance is given by

What can we say about its distribution?

Albanese and Lawi (2005) consider the Laplace transform of such integrals (see also Hurd and Kuznetsov (2008)) in general, and in particular for the CIR/Heston-process:

where and are two Borel functions. However, they come to the conclusion in Corollary 3, eq. (50), that the Laplace transform of the integral in our case, which includes an exponential factor, is not computable in closed form.

Next, we examine the -factor model, in which futures returns follow the SDE

(22) |

and the two variance processes follow the SDEs

(23) | ||||

(24) |

The correlations are given by , , and all other correlations are zero.

Inserting (22) into (16) gives for the instantaneous covariance

(25) |

In contrast to the -factor model, the instantaneous correlation in the -factor model is now stochastic. The same holds of course for the general multi-factor model with .

What can we say about the distribution of ? It follows from the definition (17) that , so that the returns of the two futures contracts are always positively correlated. To get some more insight, we consider the -factor model in a numerical example. The parameters of the model have been chosen for illustrative purposes and are given in Table 1: the first factor is more volatile than the second one, and it also decays more slowly.

For two contracts with maturities and years, respectively, we plot the empirical density function of in Figure 1.

These plotted empirical probabilities were obtained by sampling (17) one million times in a Monte Carlo simulation. The empirical mean is . In case both stochastic volatilities are made deterministic by setting in equations (23) and (24), the empirical mean is , which is in excellent agreement with the deterministic instantaneous correlation of a corresponding -factor Clewlow-Strickland model with volatility functions with

## 3 Calendar Spread Options and Analysis of Dependence

In this Section we review the definition, functioning and pricing of calendar spread options. We then review the notion of implied correlation associated to a price of calendar spread option. Finally we introduce analytic results to obtain the copula function and the copula density produced by a model defined by means of its joint characteristic function.

### 3.1 Calendar Spread Options written on WTI futures

Calendar spread options (CSO) are very popular options in commodities markets. There are two types of these options: calendar spread calls (CSC) and calendar spread puts (CSP). Like spread options in equities derivatives markets, their payoff depends on the price difference of two underlying assets. A call spread option on two equity shares and gives the holder, at time , the payoff and a put the payoff In the case of calendar spread options, the two underlyings are two futures contracts on the same commodity, but with different maturities and . Along with the volatilities, the dependence between the two contracts has a large influence on the option’s price. Note that CSOs written on commodities futures should not be confused with the well-known options strategy named calendar spread. This strategy involves two vanilla options (one bought and one sold) with different maturities, whereas the CSO is a single option. Examples of CSOs are the NYMEX calendar spread options on crude oil (WTI). A WTI CSC (CSP) represents an option to assume a long (short) position in the first expiring futures contract in the spread and a short (long) position in the second contract. There are also so-called financial CSOs traded on the NYMEX, which are cash settled. For pricing purposes we will not distinguish between these two settlement types in this paper.

There is usually very good liquidity on -month spreads (for which month), whereas options on and -month spreads are less liquid. The NYMEX CSO on -month WTI spreads can be accessed in Bloomberg using the ticker WA.

Let two futures maturities , an option maturity , and a strike (which is allowed to be negative) be fixed. Then the payoffs of calendar spread call and put options, and , are respectively given by

(26) |

(27) |

To evaluate such options with a pricing model, the discounted expectation of the payoff must be calculated in the risk-neutral measure. Assuming a continuously-compounded risk-free interest rate , we have at time :

(28) | ||||

(29) |

Note that there is a model-independent put-call parity for calendar spread options:

(30) |

Apart from Monte-Carlo simulation (where simulation of the CIR/Heston process is well-understood), we are aware of three efficient methods to price spread options. The first two are suitable when the joint characteristic function is available. The third one is more direct but needs the marginals and joint distribution function of the underlying futures.

The formula of Bjerksund and Stensland (2011) for a joint Black-Scholes model is generalized by Caldana and Fusai (2013) to models for which the joint characteristic function is known. Strictly speaking, these methods give a lower bound for the spread option price. However, our tests lead us to agree with the above authors that this lower bound is very close to the actual price (typically the first three digits after the comma are the same), and we therefore regard this lower bound as the spread option’s price itself. Furthermore, in case the formula is exact (exchange option case).

Let be the joint characteristic function of the logarithms of two futures prices as given in equation (7). Following Caldana and Fusai (2013), the price of the calendar spread option call with maturity and strike is given in terms of a Fourier inversion formula as

(31) |

where

and

The parameter controls an exponential decay term as in Carr and Madan (1999). We found a value of to perform well in numerical applications. This method appears to be the most suitable to our model and setup.

An alternative method that also works with the joint characteristic function of the log-returns has been proposed by Hurd and Zhou (2010). In their paper, the transform of the calendar spread payoff function with a strike of is analytically calculated. The price of the corresponding option is then deduced from this analytical result. Let denote the time futures log-price. Following Hurd and Zhou (2010), the calendar spread call option with maturity , underlying futures maturities , and strike , can be priced as:

(32) |

where

and is the complex gamma function defined for by the integral . The double integral in (32) is evaluated numerically using the two-dimensional Fast Fourier Transform ( FFT). The algorithm returns a whole matrix of option prices at different values of and . Options with other strikes () are then evaluated by re-scaling and interpolation, if necessary, using the same matrix.

Methods working with distribution functions instead of characteristic functions are also available to price calendar spread options. In this category of methods, the most direct approach is to evaluate a double integral of the payoff function times the joint density of the two underlying futures contracts. However, following Tavin (2014), we can write calendar spread option prices as single integrals over the marginal and joint distribution functions. The calendar spread call and put option prices are given, at and for , by

(33) | ||||

(34) |

where and are the marginal distribution functions of and , respectively, and is their joint distribution function. The case is treated as a calendar spread option written on the reverse spread with the opposite strike . A proof for spread options that can easily be adapted to our case is given in Tavin (2014).

In our model, the distribution functions involved in (33) and (34) are not readily available. However, it is possible to calculate and from the joint characteristic function of () using direct inversion formulas given by

(35) | ||||

(36) |

with a proper choice of the smoothing parameter . We find that works well in our applications. A detailed proof of these inversion results can be found in Courtois and Walter (2014). The joint distribution function can be recovered in a similar way using a direct two-dimensional inversion formula.

###### Lemma 3.1

(37) |

### 3.2 Implied Correlation for Calendar Spread Options

Given the observed price of a calendar spread option it is possible to extract an implied quantity reflecting the implied level of dependence embedded in the given price. This quantity is named implied correlation and can be defined as the parameter of the Gaussian copula that reproduces the observed price. This definition has the advantage to be well defined in the sense that implied correlation exists as soon as the observed price is free of arbitrage. This copula-based definition also has the advantage to disentangle the impact of marginals from the dependence structure on the price of a calendar spread option. This notion of implied correlation is introduced and detailed in Tavin (2014).

For , we denote by the bivariate Gaussian copula with parameter . With special cases, for , and , where and are the usual upper and lower Fréchet-Hoeffding bounds. For definitions and general theory about copula functions we refer to Nelsen (2006) and Mai and Scherer (2012).

When the chosen dependence structure is given by a Gaussian copula with correlation parameter , and its marginals are and , the price of the strike calendar spread call option is denoted by and is given by, for ,

(38) |

and, for

where and denote the prices obtained for the calendar spread option when the chosen dependence structures are respectively and .

When a calendar spread option quote is observed, it is possible to extract the implied correlation as the value of the correlation parameter to be used in (38) such that it reproduces the given price. The retained notion of implied correlation does not rely on a bivariate extension of the Black-Scholes-Merton model. Instead, it corresponds to a bivariate Gaussian copula that pairs the true underlying marginals so that the implied correlation does not imbed errors made on the marginals when using a joint log-normal model.

Let be an observed arbitrage-free price for the calendar spread option call written on and , with maturity and strike . The associated implied correlation exists and is unique. It can be obtained by solving numerically the equation in that is written

(39) |

For proofs and more details on the computation of implied correlation we refer to Tavin (2014). Note that, here, the notion of implied correlation has been reviewed with calendar spread calls. The definition can be given identically for a put option as, by put-call parity, the value is the same as for the calendar spread call option with same characteristics. On the WTI calendar spread options market, implied correlation depends on both the strike and the maturity of options. By analogy with implied volatility, these phenomena are referred to as implied correlation smile (or frown) and implied correlation term-structure.

### 3.3 Analysis of the Dependence Structure Between two Futures

As for distribution functions in our model, the dependence structure between two futures with given maturities is not readily available. In order to analyze the dependence between futures prices at a future time horizon created by our model we need to work with the joint characteristic function. As we have seen, it is possible to recover marginal and joint distribution functions directly from the joint characteristic function. It is also possible to recover via analytical expressions the copula function and the copula density that characterize the dependence between two futures at the chosen time horizon. These analytical expressions are also valid in the more general context of a two-dimensional stochastic process defined by means of its joint characteristic function.

Let with the chosen time horizon and , the expiry dates of the pair of futures. The joint density