Calendar Spread Exchange Options Pricing with Gaussian Random Fields

: Most of the models leading to an analytical expression for option prices are based on the assumption that underlying asset returns evolve according to a Brownian motion with drift. For some asset classes like commodities, a Brownian model does not ﬁt empirical covariance and autocorrelation structures. This failure to replicate the covariance introduces a bias in the valuation of calendar spread exchange options. As the payoff of these options depends on two asset values at different times, particular care must be taken for the modeling of covariance and autocorrelation. This article proposes a simple alternative model for asset prices with sub-exponential, exponential and hyper-exponential autocovariance structures. In the proposed approach, price processes are seen as conditional Gaussian ﬁelds indexed by the time. In general, this process is not a semi-martingale, and therefore, we cannot rely on stochastic differential calculus to evaluate options. However, option prices are still calculable by the technique of the change of numeraire. A numerical illustration conﬁrms the important inﬂuence of the covariance structure in the valuation of calendar spread exchange options for Brent against WTI crude oil and for gold against silver.


Introduction
Gaussian fields have been used for several decades in the analysis of spatial statistics. One of the key features in spatial statistics is the autocorrelation of data. The observations at locations in close spatial proximity often tend to be more similar than observations at locations far apart. We refer the interested reader to the books of Cressie (1993), Adler (1981) or Matern (1986) for a detailed presentation of spatial statistics and the theory of Gaussian fields. An introduction to Gaussian processes on general parameter spaces may be found in Adler and Taylor (2009).
We also find applications of Gaussian fields in the financial literature, e.g., Goldstein (2000) generalized the work of Kennedy (1994) and proposed a model for interest rates based on a two-dimensional random field. In this model, increments along time are independent, but the correlation structure between bond yields of different maturities can be arbitrarily chosen. The work in Kimmel (2004) introduced a state-dependent volatility in this model. Albeverio et al. (2004) extended previous models for yield curves with Lévy fields. The work in Ozkan and Schmidt (2009) used similar fields to model the yield curve in the presence of credit risk.
Gaussian fields are also used for actuarial modeling. The work in Biffis and Millossovich (2006) modeled the intensity of mortality as a random field, in order to capture cross-generation (risk class) effects induced by the on-going management of portfolios of policies. Based on this framework, Biagini et al. (2017) studied the pricing and hedging of life insurance liabilities. The work in Wu (2016) developed Gaussian process regression methods for mortality forecasting.
The pricing of exotic options with a payoff involving asset prices at different times requires a model able to replicate the covariance of underlying securities. Assuming that asset returns are ruled by a Brownian motion with drift is convenient for mathematical developments. However, this model does not replicate the time dependence observed for some asset classes as underlined by Willinger et al. (1999). On the other hand, option prices usually do not admit any closed form expressions for models with stronger econometric foundations. These points motivate this research that has a two-fold objective. Firstly, we propose a Gaussian model that is analytically tractable, with three possible covariance structure for asset prices. Secondly, we aim to emphasize the importance of the covariance in the valuation of exotic derivatives. For this purpose, we derive closed form expressions for calendar exchange spread and Asian calendar exchange spread options. Next, we evaluate in numerical applications the bias induced by a misspecification of the covariance structure.
In our approach, asset prices are seen as Gaussian fields indexed by time. These fields are conditional in the sense that their initial value is known. We work with homogeneous fields: the covariance only depends on the length of the time interval between considered prices. We study three covariance structures with sub-exponential, exponential and hyper-exponential autocorrelation functions. After a presentation of model features, we propose a method to simulate conditional Gaussian fields based on the spectral decomposition of the covariance function. Within this approach, a simulated sample path is known continuously at all times. This characteristic is particularly useful for pricing path-dependent options.
In order to emphasize the role played by the correlation on option prices, we evaluate two exotic derivatives with a payoff depending on asset prices at different times. The first one is the calendar spread exchange option. This option pays the positive difference between two asset prices observed at different times. The second derivative is an Asian calendar spread exchange option with geometric averages. This product pays the positive difference between the geometric average returns of two assets, calculated over different time intervals. Asian options with arithmetic averages are not studied as their price does not admit any closed form expression. In the numerical illustration, we calibrate our models to two pairs of assets (Brent-West Texas Intermediate (WTI) crude oil, silver-gold) and evaluate exchange options. Our results confirm that the misspecification of the autocorrelation function has a strong impact on the prices of calendar spread exchange options.

Conditional Gaussian Fields
We consider a probability space (Ω, F , Q) endowed with the filtration (F t ) t≥0 on which is defined a random field X t indexed by time. Before detailing the features of asset prices, we recall in this section the main properties of Gaussian random fields. A real valued random field is a family of random variables X t indexed by t ∈ R + with a collection of distribution functions of the form F t 1 ,...,t n : F t 1 ,...,t n (b 1 , . . . , b n ) = P (X t 1 ≤ b 1 , . . . , X t n ≤ b n ) , where b 1 , ..., b n ∈ R. In this article, we consider Gaussian random fields for which the distribution F t 1 ,...,t n (.) is a multivariate normal distribution with a zero mean, E [X t ] = 0. This distribution is fully characterized by its covariance matrix. The covariance between X t and X s is denoted by C(t, s). This is a non-negative definite function on (R + × R + ). Given that the mean of the Gaussian field is null, the covariance function is equal to its cross expectation: If the covariance function is equal to C(t, s) = t ∧ s, the Gaussian field is a pure Brownian motion with independent increments. However, a Gaussian field does not have in general independent increments and is not a Markov process. On the other hand, a Gaussian field is not systematically a semi-martingale. If the Gaussian field is not a semi-martingale, we cannot rely on stochastic differential calculus for evaluating options. In counterpart, working with Gaussian fields allows us to reproduce with accuracy the autocovariance function of asset prices compared to processes with independent increments. We will see later that this feature is particularly important for pricing calendar options. In this article, we consider continuous and differentiable covariance functions, C(t, s) ∈ C ∞ , and focus on the homogeneous Gaussian random field. A random field is homogeneous or stationary if E X 2 t is finite for all t and: 1.
E [X t ] is constant and independent of t ∈ R + , 2.
C (t, s) solely depends on the difference |s − t|.
When the Gaussian field is homogeneous, the covariance function may be rewritten as a function g : R → R + , The construction of the function g(.) is detailed in the next section. For homogeneous fields, increments in the t direction are correlated. For t, h ∈ R + , the time covariance is equal to: This is a clear difference from the Brownian motion, which has independent increments. Furthermore, the conditional expectation of X s with respect to the filtration F t depends on the whole sample path of the process and not exclusively on X t . In general, X t is then not a Markov process 1 .
In financial applications, X 0 is known at Time 0 and will determine the asset prices. Therefore, we must determine the conditional distribution of X t with respect to X 0 . Given that (X t , X 0 ) is a bivariate Gaussian distribution: The random variable Z t := X t |X 0 = x is Gaussian d ∼ N µ t|x , σ 2 t|x with a mean and variance given by: Using the properties of the multivariate normal distribution, we directly infer the covariance between Z t and Z s in the next proposition that is used in later developments. Proposition 1. The joint distribution of Z t = X t |X 0 = x and Z s = X s |X 0 = x for 0 ≤ t ≤ s is a bivariate normal random variable with a mean and covariance matrix given by: Excepting the Ornstein Uhlenbeck process, that is a Markov process that may be reformulated as an homogeneous Gaussian field. However, in general, a time Gaussian field is not a Markov process.
In the next section, we consider covariance functions of homogeneous fields that are decreasing functions of h such that lim h→∞ g(h) = 0. From the previous proposition, we infer the asymptotic variance and covariance of Z t : and: The variance being bounded by a constant, the process Z t reverts around its mean. If X t is a Brownian field, the covariance between X t and X s for 0 ≤ t ≤ s is equal to C [X t , X s ] = s ∧ t = t. Conditional to X 0 , Z t = X t |X 0 = x and Z s = X s |X 0 = x are distributed according to a bivariate normal distribution with a mean (x, x) and a covariance matrix s t t t . In this particular case, when t → ∞, the variance tends to infinity, and the covariance is equal to t: Therefore, homogeneous fields have a very different behavior from the Brownian motion. As their variance is asymptotically constant, homogeneous fields revert to their mean and are not adapted to model stock prices, which usually have a variance increasing with the time horizon. This property makes them more suitable for the modeling of interest rates or any assets with a mean reverting yield. The proposed approach is also appropriate for commodities; mainly because their value reverts to a long run equilibrium level with a constant asymptotic variance as illustrated in Bessembinder et al. (1995). Gibson and Schwartz (1990); Cortazar and Schwartz (1994); Schwartz (1997) or Hainaut (2017) modeled commodities with an Ornstein-Uhlenbeck (OU) process presenting the same features as those of Equation (1).

Market Model
The filtration (F t ) t≥0 is generated by d independent continuous and homogeneous random fields X (k) t for k = 1 to d. They have a null mean, and the covariance function of X (k) t is denoted by g k (h) := C k (t, t + h) where g k : R → R + and g k (h) = g k (−h). The processes X (k) t are homogeneous in time, but increments of the random field in the t direction are correlated: for t, h ∈ R + , we have that C X The probability measure Q is here the probability measure under which the pricing of financial derivatives is done.
As in the previous section, we denote by Z The risk-free rate is assumed constant and is denoted r. We assume that the financial market counts two securities, S (1) t and S (2) t , with a homogeneous covariance structure. The yields 2 of these assets are denoted by y (1) t and y (2) t . The relation between yields and prices is S The yield y (j) t is here defined by the cumulative return up to time t of the j th asset. j = 1, 2. We denote the vector of conditional Gaussian fields by . In our model, the yields are linear functions of Gaussian fields: where α j (t) for j = 1, 2 are functions from R + to R + equal to: Σ is a 2 × d real matrix 3 . c 1 and c 2 are adjustments of the drift for non-financial assets like commodities. Under the risk-neutral measure, commodities earn on average the risk-free rate plus the cost of carry and less the convenience yield. If S (1) t and S (2) t are financial assets, these factors are null: c 1 = c 2 = 0. According to the usual financial theory, discounted financial security prices must be martingale processes under the pricing measure Q. More precisely if Σ j,. is the j th line of Σ, discounted prices are such that E e −rT S (j) 0 E e −α j (T)+Σ j,. Z T |F t for any times t ≤ T. Except in the case of the Brownian motion, a time Gaussian field depends in general on its sample path in a complex manner, and the process is not Markovian. Therefore, the expectation E e Σ j,. Z T |F t does not admit a closed form solution, and the drift α j (T) that ensures that discounted prices are always martingale is unknown 4 . However, at time t = 0, the martingale condition E e −rT S (j) (2) t are non-financial assets, their future price is equal to E S (j) 0 e (r+c j )T for j = 1, 2, and this relation holds only if α j (t) are defined by Equation (3). If we remember Proposition 1 and given that g k (0) = 1, the autocovariance of y (j) t is: whereas the covariance of y (1) t and y (2) s is equal to: We remark that this is an important difference with existing models using Gaussian fields, e.g., Goldstein (2000) or Kennedy (1994) developed a model for interest rates in which time increments are independent. In this particular case, the Gaussian field is a semi-martingale, and we can use the Itô's calculus to deduce the main properties of asset dynamics. 3 WE remark that Σ may eventually be replaced by a time-dependent matrix Σ(t) in order to replicate seasonality effects in the covariance. In this framework, the pricing by simulations would still be possible.

Choice of Autocovariance Functions
In this article, we consider three types of autocovariance functions that are built with Bochner's theorem. A function g(h) is eligible to be an autocovariance function if it is positive definite. The work in Bochner (1993) showed that positive definite functions may be defined by their spectral measure.
Theorem 1. (Bochner) A continuous function g(h) from R + to the complex plane C is positive definite if and only if it may be represented as the Fourier transform of a measure ν(.) on R: where ν(.) is a bounded, real valued function such that A ν(du) ≥ 0 for all A ⊂ R.
The function ν(.) is called the spectral distribution function for g. An alternative representation is: Then, if Z t is R-valued, g(h) is also R-valued, and the imaginary part in this last equation is null: R sin (hu) ν(du) = 0. This implies that the spectral distribution is symmetric around the origin. If ν(.) is continuous, we then have ν(du) = ν(−du) and: In this work, we consider three spectral distributions ν(.): a double exponential, an inverse quadratic and a double quadratic exponential measure. The next paragraphs detail each one.
(1) The first measure considered is the double exponential measure. Let λ ∈ R + be constant. A symmetric exponential measure is defined by the following relation: If we use this function as the spectral distribution of the covariance, we obtain the following function g(h): which is a sub-exponential function (in the sense that it decays at a lower pace than an exponential function). In order to define autocovariance functions of Gaussian fields Z (k) t involved in the dynamics of the financial market, we consider d spectral density functions (4) parameterized by a vector λ = {λ 1 , . . . , λ d } ∈ R d+ . The corresponding autocovariance functions are in this case denoted g sub k (.) and defined by: (2) The second type of spectral measure that we study is the inverse quadratic one: where λ ∈ R + . By construction, this measure is symmetric, and the function g(h) is an exponential decreasing function: It may be shown that this autocovariance function is closely related to the one of an Ornstein-Uhlenbeck process. Details are provided in Appendix A. In later developments, we consider d exponential autocovariance functions for Z (k) t parameterized by λ = {λ 1 , . . . , λ d } ∈ R d,+ and denoted by: (3) The last category of covariance functions is generated with what we call an exponential quadratic spectral measure: Given that R + e −λu 2 cos (hu) du = 1 2 π λ e − h 2 4λ , the autocovariance function is in this case equal to: We abusively call this covariance function the exponential quadratic covariance. This function decays at a faster pace than an exponential one and is in this sense hyper-exponential decreasing.
In the following sections, we consider d exponential quadratic covariance functions for Z (k) t parameterized by λ = {λ 1 , . . . , λ d } ∈ R d,+ and denoted as: We will price calendar spread exchange options using these three autocovariance functions g sub k (h), g exp k (h) and g qua k (h). As autocovariance functions g sub k (h) and g qua k (h) decrease respectively at a lower and at a faster pace than g exp k (h), our approach covers a wide spectrum of autocorrelation structures. Options are priced under the risk-neutral measure, and autocovariance functions should therefore duplicate option market prices. However, if the option market is not liquid enough, an alternative solution to calibrate g k (.) consists of assuming that covariances under the pricing and real measures are similar (eventually adjusted by a risk premium). Under this assumption, the functions g k (.) ∈ {g sub k (.), g exp k (.), g qua k (.)} can easily be calibrated with the method of moments matching. In this case, we minimize the quadratic spread between marginal and empirical covariances denoted bŷ C [X . , X .+∆h ]: where ∆ is the time interval between two successive observations. This approach is applied in numerical illustrations in order to determine which covariance function is the most suitable for some commodity prices (silver, gold, Brent and WTI crude oil).

Simulation of a Conditional Field by Spectral Decomposition
In this section, we propose a method to simulate sample paths of asset prices. As increments of Gaussian fields involved in the dynamics of S (1) t and S (2) t are not independent, the simulation of these processes requires particular care. A natural approach consists of simulating the asset dynamics with a multivariate normal distribution conditioned by X 0 = 0, on a discrete time grid. However, there exists an elegant alternative based on a discretization in the space of frequencies. With this approach, the sample path of S (1) t and S (2) t is known at all times t and not only at discrete times. This feature is particularly interesting for the pricing of path-dependent options. Let us first recall that the covariance function of a homogeneous Gaussian field is defined by a measure ν(.). In the developments of Section 4, we have considered a (probability) measure 5 ν(.) of R. For any Borel subset A ⊂ R, we note ν(A) = A ν(du). We define a complex noise W based on the measure ν, or ν-noise, as a random process defined on Borel subsets of R, such that for all A, B ∈ B with ν(A) and ν(B) finite, we have: When the measure ν(.) is the Lebesgue measure and W(A) d ∼ N (0, ν(A)), W is a white noise and more precisely a Brownian motion. Having defined complex ν-noises, we build the integral with respect to W, of a function ϕ(t) : R → C such that ϕ(x) 2 ν(dx) < ∞. By analogy with the Riemann integral, the integral of ϕ(t) over an interval R with respect to W is defined as the limit of a sum over a By construction, R ϕ(u)W(du) is a random variable that has zero mean and variance given by This construction allows us to state a useful representation theorem: 5 Notice that ν(.) can be any measure on R, of probability or not (e.g., Lebesgue measure). Theorem 2. (Spectral representation theorem) Let ν be a finite measure on R and W complex ν-noises. Then, the complex valued random field: has a covariance function: If W is Gaussian, then so is X t . Furthermore, to every mean-square centered (Gaussian) stationary random field on R, with covariance function C and spectral measure, there corresponds a complex (Gaussian) ν−noise W on R such that Equation (10) holds in the mean square for each t ∈ R N . In both cases, W is called the spectral process corresponding to X t .
In one direction, the proof of this theorem is immediate. It is a consequence of the construction of the stochastic integral W(ϕ) that X t defined by Equation (10) has covariance function (11). The other direction is less direct, and we refer to the book of Adler and Taylor (2007) for a proof (Theorem 2).
In this article, we consider a real Gaussian field X t with a symmetric spectral measure. Then, its covariance function can be rewritten as follows: where we have defined a positive measure µ(A) = 2ν(A) for all A ⊂ R + . In Section 4, we have introduced three symmetric measures leading to three types of covariance. The measure µ(.) for each of these cases is given by: The next proposition shows that we can reformulate the Gaussian field as integrals with respect to µ-noises: Proposition 2. A Gaussian field X t with a covariance function (12) can be rewritten as the sum of two integrals with respect to independent real (Gaussian) valued µ−noises W 1 and W 2 such that: The proof is reported in Appendix B. In practice, we use Equation (13) to approximate the sample t over a partition {u 0 , u 1 , ..., u n } of R + (with u 0 = 0) as follows: where W (k) 1 (.) and W (k) 2 (.) are normal random variables: If the measure is double-exponential (and the function g is sub-exponential), then the measure of a time interval is equal to: If the measure is inverse quadratic (and the function g is exponential), then the measure of a time interval is equal to: If the function g is exponential quadratic, then: The number of intervals n and u n is chosen such that the variance 2 , is close to zero. Equation (14) allows us to simulate a Gaussian field X t that is not conditioned by its initial value X 0 . In order to simulate Z (k) and only if the following random variable is null: By construction, the random vectorW where 0 n+1 is a null vector of dimension n + 1 and Σ k is a (n + 1) × (n + 1) matrix of covariance: 1 , H is a multivariate normal distribution with the following mean and covariance: Using the properties of the conditional Gaussian multivariate distribution, we have thatB (k) := with a covariance matrix: and where 1 n+1 Σ k Σ k 1 n+1 is a n + 1 × n + 1 matrix: In numerical illustrations, the sample path of Z (k) t over a partition {u 0 , u 1 , ..., u n } of R + (with u 0 = 0) is then simulated with the following relation: and W (k) 2 (.) are normal random variables: This technique allows for simulating the sample path of Z t , which are known for all times t ∈ R + . This method is particularly helpful for the pricing of path-dependent options.

Calendar Spread Exchange Options Pricing
One objective of this article is to emphasize the importance of the autocorrelation in the valuation of exotic derivatives. For this purpose, we derive closed form expressions for calendar spread exchange options and for Asian options of the same type in the next section. Next, we will assess numerically the bias induced by a misspecified covariance structure. A calendar spread exchange European option is a financial derivative delivering a payoff equal to the positive difference between β 1 S (1) T and β 2 S (2) t , at expiry T. The price of this option is equal to the expected discounted payoff under the risk-neutral measure: Gaussian fields Z (k) t have sub-exponential, exponential or exponential quadratic autocovariance functions: g k ∈ {g sub k (.), g exp k (.), g qua k (.)} for k = 1, ..., d. In order to obtain a closed form expression for the option price, we introduce a new probability measure, denoted by Q S 2 . This measure uses S (2) t as numeraire. The change of measure from Q to Q S 2 is defined by the following random variable: where α 2 (t) = 1 2 ∑ d k=1 (σ 2k ) 2 1 − g k (t) 2 and Σ 2,. is the second line of the matrix Σ. This approach was used by Margrabe (1978) so as to evaluate the exchange option in a Brownian setting. By definition, dQ S 2 dQ is a strictly positive, F t -measurable random variable such that E dQ S 2 dQ |F 0 = 1. Conditional to F 0 , it defines then an equivalent measure Q S 2 to Q. For any F T -adapted random variable N T , we have the following relation: We immediately infer from this last expression that the expectation of N T under Q is also equal to: This result allows us to rewrite the calendar spread exchange option as follows: (17) The next step consists of determining the statistical distribution of By construction, the ratio of asset prices is given by: where Σ 1,. and Σ 2,. point out the first and second line of the matrix Σ. Let us denote the logarithm of this ratio by: The statistical distribution of R t,T under the measure Q S 2 is Gaussian, as detailed in the next proposition proven in Appendix B : Proposition 3. R t,T under the measure Q S 2 is a normal random variable with a mean µ R (t, T) and a variance σ R (t, T) 2 equal to: The exchange option price as reformulated in Equation (17) is a call option on a lognormal underlying asset. As stated in the next proposition, the option price admits then a closed form expression similar to the Black and Scholes equation, and this result holds whichever the covariance function g k ∈ {g sub k (.), g exp k (.), g qua k (.)}.
Proposition 4. Let Φ (x) be the cumulative distribution function of a standard normal random variable. The price of a calendar spread exchange option price is given by the following expression: where d 1 (t, T) and d 2 (t, T) are defined as follows: and: In Section 8, we evaluate calendar spread exchange options for Brent against WTI crude oil and for silver against gold.

Asian Calendar Spread Exchange Options, with the Geometric Average
Our approach presents the same level of tractability as a Brownian motion. To illustrate this point and to emphasize the role of the covariance structure on Asian options, we price another type of exotic derivative that depends on the whole sample path of asset prices. We focus on Asian calendar spread options with a payoff related to geometric average returns of assets S T and β 2 G (2) t , at expiry T. The price of this option is the expected discounted payoff under the risk-neutral measure: By definition, the log-return of assets is equal to: Therefore, the integral of this log-return over the interval [0, t], which is equal to ln G Before any further developments, we present the properties of the integral of Z Given that Z (k) t j is a Gaussian random variable, this integral is distributed as a normal random variable of a null mean under the condition that its variance exists. The variance of T 0 Z (k) u du is obtained by considering the limit of the product of integrals over the partition Π n : t for k = 1, ..., d are independent processes, the expectation of this product with respect to F 0 is equal to the variance of This last double integral is bounded by T 2 g(0) for the three covariance functions that are considered in this article, and the integral From the definition (21) From Equation (22) and Proposition 5, we immediately infer the next corollary: Proposition 6. If covariance functions are sub-exponential, g k (h) = g sub k (h), the variance of ln G (j) t , denoted by σ ln G j (t) 2 , is equal to: The drift of ln G (j) t depends on the integral of α j (.). The next proposition proven in Appendix B provides the analytical integral of this function when the covariance is sub-exponential.
Proposition 7. If covariance functions are sub-exponential, g k (h) = g sub k (h), then the integral of function α j (.) is given by: If the covariance function is exponential, then the double integrals needed to evaluate the variance of t 0 Z s ds in Equation (22) are given by the following proposition detailed in Appendix B.
whereas the double integral t 0 T 0 g k (u)g k (v)dvdu is given by: A direct corollary of this proposition is that the variances of ln G (j) t for j = 1, 2 admit a closed form expression: Proposition 9. If autocovariance functions are exponential g k (h) = g exp (h), the variance of ln G (j) t , denoted by σ ln G j (t) 2 , is equal to: for j = 1, 2.
The drift of ln G (j) t depends on the integral of α j (.) which is provided in the next proposition, proven in Appendix B.
When covariance functions are exponential quadratic, the integrals involved in the variance of ln G (j) t are given in the following proposition: Proposition 11. Let ω ∈ C − and erf(x) = 2 √ π x 0 e −t 2 dt be the error function. If g k (h) = g quad k (h), then: When the covariance functions are quadratic exponential, the variance of ln G (j) t for j = 1, 2 is given by the next proposition: Proposition 12. If the autocovariance functions g k (h) are defined by Equation (8), then the variance of ln G (j) t , denoted by σ ln G j (t) 2 is equal to: for j = 1, 2.
The drift of ln G (j) t depends on the integral of α j (.). The next proposition proven in Appendix B provides the analytical integral of this function when the covariance is quadratic exponential.
In order to price the Asian calendar spread exchange option, we introduce a new probability measure, denoted by Q G 2 . The change of the measure from Q to Q G 2 is defined by the following random variable: By definition, dQ G2 dQ is a strictly positive random variable, F t -adapted and such that E dQ G 2 dQ |F 0 = 1. It defines then an equivalent measure Q G 2 to Q. For any F T -adapted random variable N T , the following relation holds: The expectation of N T under the risk-neutral measure is then equal to: t is a lognormal random variable. Using this change of measure allows us to develop the option price as follows: Next, we study the moment generating function of the ratio . If we remember the definitions of G (1) T and G (2) t , this ratio may be developed as follows: where the exponent in this last expression is a random variable denoted by A t,T and defined as follows: The next result proven in Appendix B states that A t,T is a normal random variable under Q G 2 . (33) is a normal random variable with a mean and variance respectively equal to:

Proposition 14. The random variable A t,T defined in Equation
and: Finally, we infer a closed form expression for the Asian exchange calendar option. The proof is detailed in Appendix B.
Proposition 15. Let Φ (x) be the cumulative distribution function of a standard normal random variable. The value of an Asian exchange calendar option is given by the following expression: where: and: In the next section, we compare Asian option values obtained with the different covariance structures fitted to some metals and oil prices.

Numerical Illustration
The market of exotic options lacks liquidity, and the absence of data prevents us from estimating models by replicating market prices. An alternative solution for calibrating functions g k (.) consists of assuming that the covariance structures under the pricing and real measures are similar (eventually adjusted by a risk premium). We adopt this approach and calibrate models by minimizing the quadratic spread between marginal and empirical covariances. The model is fitted to two pairs of commodities 6 : Brent vs. WTI oil and silver vs. gold. The dataset contains daily log-returns from 11 April 2008-11 April 2018 (2522 observations). The number of lags considered in the optimized objective (9) is set to 150 days. The square roots of the quadratic error between empirical and model covariances are reported in Tables 1 and 2 for different numbers of homogeneous fields. We draw two conclusions from these figures. Firstly, increasing the number of fields reduces the error whatever the covariance 6 As mentioned in Section 2, homogeneous fields have constant asymptotic variance. This feature makes them more suitable for the modeling of commodities or interest rates. structure. Secondly, the best fit is obtained with sub-exponential covariance functions. Tables 3 and 4 contain parameter estimates with d = 3 fields.  Recall that the field with an exponential covariance corresponds to an Ornstein-Uhlenbeck (OU) process. Modeling commodity prices with this mean reverting process as in Cortazar and Schwartz (1994) is mathematically convenient; mainly because we can rely on stochastic calculus to price options. However, our results clearly suggest that these processes do not replicate the empirical covariance as well as sub-exponential Gaussian fields. This is confirmed by Figure 1, which compares the covariance structures of fitted models with d = 3 homogeneous fields to empirical covariances. The observed autocovariances are concave decreasing functions, and the exponential model fails to replicate this concavity, particularly in the short term.  Tables 5 and 6 report calendar spread option prices 7 written on the pairs Brent-WTI crude oil and silver-gold. These figures emphasize the impact of a misspecified autocovariance function on the value of calendar spread options. Assuming an OU dynamics for Brent and WTI oil leads to overestimating option values by 5-37% compared to those computed with a sub-exponential covariance. The gap between option prices computed with quadratic exponential and sub-exponential covariances is much smaller and ranges from −1.55-2.95%. For silver against gold, the conclusions are similar: the bias caused by the wrong choice of covariance structure can be significant (up to 6.2% with an OU model and −19.6% with a quadratic exponential covariance). Table 7 presents calendar spread option prices computed by Monte Carlo simulations with different levels of discretization in the space of frequencies.
As expected, increasing u n and decreasing the size of the step of discretization reduce the gap between values obtained analytically, as well as the numerical estimates.

Calendar Spread Exchange Options
t ) + where T = 1 year and t ranges from 0.3-0.9 years. β 1 and β 2 are such that β 1 S (1) To understand the influence of autocorrelation parameters λ k on option values, we perform a sensitivity analysis for the sub-exponential model driven by d = 2 conditional Gaussian fields. Table 8 reports the parameters of this model fitted to Brent and WTI crude oil. We appraise calendar spread exchange options with three sets of autocovariance parameters: (λ 1 , λ 2 ), λ 1 2 , λ 2 and λ 1 , λ 2 2 . Figure 2 emphasizes the impact of these modifications on autocovariance functions g 1 (t) and g 2 (t). We observe that reducing λ 1 or λ 2 decreases the level of autocorrelation of Z (1) t and Z

Calendar Spread Exchange Options
(2) t . Table 9 underlines the important role played by λ 1 or λ 2 on option prices. The relative spread between prices computed with (λ 1 , λ 2 ) and λ 1 2 , λ 2 ranges from 5% up to 55%; whereas, the spread between prices 7 The discount rate is set to r = 5% and c 1 = c 2 = 0. obtained with (λ 1 , λ 2 ) and λ 1 , λ 2 2 varies from 7.5-11%. These results confirm the importance of the autocorrelation in the valuation of exotic derivatives like calendar spread exchange options.

Sub-Exponential Parameters for the Sensitivity Analysis
(2) t ) + where T = 1 year and t ranges from 0.3-0.9 years. β 1 and β 2 are such that β 1 S (1) Tables 10 and 11 report Asian calendar spread option prices on Brent vs. WTI crude oil and silver vs. gold. The spreads between prices obtained with sub-exponential and exponential covariance structures are smaller than those of calendar options. However, they remain significant: from −1.84-3.48% for Brent vs. WTI and from −1.74-1.54%. Prices obtained with an exponential quadratic covariance model deviate more widely from those computed with a sub-exponential covariance (from 7.16-20.6%).

Asian Calendar Spread Exchange Options
t ) + where T = 1 year and t ranges from 0.3-0.9 years. β 1 = β 2 = 1. To understand the influence of autocovariance parameters λ k on Asian options, we perform a sensitivity analysis for the sub-exponential model driven by d = 2 conditional Gaussian fields. We appraise Asian calendar spread exchange options with the sets of autocovariance parameters, (λ 1 , λ 2 ), λ 1 2 , λ 2 and λ 1 , λ 2 2 , presented in Table 9. Table 12 reveals that the relative spread between prices computed with (λ 1 , λ 2 ) and λ 1 2 , λ 2 ranges from 0.75-4.3%, whereas the spread between prices obtained with (λ 1 , λ 2 ) and λ 1 , λ 2 2 varies from 0.6-5.5%. The impact of the autocorrelation on Asian option prices is still significant, but less important than for exchange options due to the smoothing induced by the definition of the payoff.

Conclusions
This article has two objectives. Firstly, it proposes a simple alternative to models based on Brownian motions that allows for an accurate modeling of covariance and autocorrelation structures. In our approach, price processes are seen as conditional homogeneous Gaussian fields indexed by time. As this field is in general not a semi-martingale, we cannot rely anymore on stochastic calculus to evaluate options. However, a conditional Gaussian field inherits all the properties of the multivariate normal random variable. This guarantees its high analytical tractability.
Homogeneous Gaussian fields indexed by time have a covariance function depending on the lag between two realizations of this field. They include Ornstein-Uhlenbeck processes, and their asymptotic variance and autocovariance are respectively constant and null. The variance of these fields being bounded by a constant, asset yields revert around their mean. This property makes the proposed approach more suitable for the modeling of interest rates or any assets with a mean reverting yield, like commodities.
We study three types of covariance structures with sub-exponential, exponential and hyperexponential decreasing autocorrelations. We also propose a procedure of simulation based on the spectral decomposition of the covariance matrix. Within this approach, the sample path of simulated prices is known in continuous time. This feature is particularly interesting for the numerical pricing of path-dependent options.
Secondly, we aim to emphasize the importance of the covariance structure in the valuation of calendar spread derivatives. We derive closed form expressions for Asian and non-Asian calendar spread exchange options. Analytical formulas are obtained by the technique of the change of numeraire. Next, we fit models to two pairs of commodities: Brent vs. WTI oil and silver vs. gold. The numerical analysis reveals that these assets exhibit concave decreasing autocovariances that the exponential model fails to replicate. For these asset classes, the best fit is obtained with sub-exponential covariance functions. We observe that covariance misspecification introduces a bias of several percent for the prices of calendar spread exchange options. The numerical illustration emphasizes that assuming an OU dynamics for Brent and WTI oil leads to overestimating option values by 5-37% compared to those computed with a sub-exponential covariance structure. For gold-silver exchange options, the bias introduced by the wrong choice of covariance structure can be significant (up to 6.2% with an OU model and −19.6% with a quadratic exponential covariance). An analysis of the sensitivity of options to a reduction by a half of the autocorrelation factors reveals that prices deviate by 5-55% from their initial values. These results confirm the importance of the autocorrelation in the valuation of calendar exchange options. For Asian options, the spreads between prices obtained with sub-exponential and exponential covariance structures are smaller than those of calendar options. However, they remain significant and range from −1.84-3.48% for Brent vs. WTI and from −1.74-1.54% for gold vs. silver. Prices obtained with an exponential quadratic covariance model deviate more widely from these computed with a sub-exponential covariance (from 7.16-20.6%). An analysis of the sensitivity of options to a reduction by a half of the autocorrelation factors reveals that option values diverge by 0.6-5.5% from their initial price. The impact of the autocorrelation on Asian option prices is still significant, but less important than for exchange options due to the smoothing induced by the definition of the payoff.
There are several interesting topics for future research. For example, it would be interesting to replace constant coefficients in Σ by a time-dependent matrix in order to model the seasonality in dependence modeling. Another improvement could be to model the risk-free rate and/or convenience yield with other Gaussian fields.
Funding: I thank the Chair "Data Analytics and Models for insurance" of BNP Paribas Cardiff, hosted by ISFA (University Claude Bernard, Lyon) and managed by the "Fondation Du Risque".
Therefore, we infer that: and X (k) 0 is distributed according to a Gaussian law, with mean zero and variance equal to one. We have that √ 2λ k t 0 e −λ k u dW (k) u is a normal random of the null mean with variance equal to: 2λ k t 0 e −2λ k u du = 1 − e −2λ k t , and e −λ k t X (k) 0 has a variance equal to e −2λ k t . The variance of X (k) t is then e −2λ k t + 1 − e −2λ k t = 1; whereas the covariance between X = 2λ k e −λ k (s+t) 1 2λ k e 2λ k (s∧t) − 1 + e −λ k (t+s) = e −λ k |s−t| .
Notice that in this article, we set X t+h = g k (s − t) − g k (s)g k (t) = e −λ k |s−t| − e −λ k (s+t) .

Appendix B. Detailed Proofs
Proof of Proposition 2. Given that W 1 and W 2 are independent and that the expectation of the integral with respect to a µ−noise is null, we have the following relation: R + sin (ut) W 2 (du) .
Proof of Proposition 4. If we rewrite the option price as an expectation under the measure Q S 2 , we obtain that: where u min = ln β 2 β 1

S
(1) 0 . According to the definition (19) of d 2 , the integral in the first term is equal to: +∞ u min f R t,T (u)du = Φ (−d 2 (t, T)) ; whereas the second integral of Equation (A4) is developed as follows: If we define a new variable v = µ R + σ 2 R , we develop this last equation as: Combining Equations (A7)-(A9) allows us to obtain Equation (26). Furthermore, we deduce by direct integration that: Proof of Proposition 11. The double integral We may check that the inner integral is given by: Furthermore, the integral of the first error function in this last equation is equal to: Combining Equations (A15)-(A16) leads to the following expression for the Asian exchange calendar option: