1. Introduction
A spread option is a financial derivative whose value depends on the difference of two underlying asset prices. When market participants are interested in a relative performance between two assets, rather than the price of one asset, the spread derivatives are suitable products for hedging or speculation. In practice, there are several types of spread options traded in the market. Certain types of spreads consist of the raw materials and the end products. A crack spread option, an option on the spread between the crude oil and the refined one, is a typical example of the production spread options. In fact, the crack spread is directly related to the profit of oil companies since it implies the cost of refining the crude oil. Among other types of production spread is the soybean crush spread. The crush means the process of converting soybean into oil and meal and so the crush spread is the difference in the values of the soybeans and the end products. The soybean crush options are listed and traded in CBOT with Globex code SOM:SI. The spark spread is the difference between the natural gas and the electricity, which reflects the efficiency of power-plant operation. Spreads can be made by futures with two different maturities on the same underlyings, which is a type of the so-called calendar spread. Many commodity products such as cotton, corn and soybeans are traded as the underlyings of such calendar spread derivatives. Also, there are spread options between similar underlyings, such as WTI-Brent oil spread futures and options. These options are listed at Inter-Continental Exchange (ICE), together with other derivatives on each of WTI oil and Brent oil.
In view of pricing spread options, the bivariate log-normal model, the two-dimensional extension of the classical Black–Scholes framework, is the basic model for practitioners. So far, there is no closed-form exact solution formula found even under this constant volatility model. The earliest research on spread options was done by Margrabe [Reference Margrabe18], which is focused on the special case with zero strike, that is, exchange options. Kirk [Reference Kirk15] suggested a closed-form approximation formula, called the Kirk formula, which is valid when the strike is non-zero but small, extending the Margrabe formula for the prices of spread options. This is the most popular formula among market participants. The Kirk formula was derived by Lo [Reference Lo17] using the idea of the WKB method. Carmona and Durrleman [Reference Carmona and Durrleman5] studied an accurate lower-bound estimation of the spread option prices. In Bjerksund and Stensland [Reference Bjerksund and Stensland2], the implicit strategy in Kirk's formula was used to obtain another approximation formula for spread options. This approximation formula is written in a closed-form expression and it is more accurate than Kirk's formula. It allows a more wide range of strike prices extending the Kirk formula. We call it the Kirk–Bjerksund–Sternsland (KBS in brief) formula in this paper. There are some recent studies on spread options. For example, Li and Wang [Reference Li and Wang16] studied spread options with counterparty risk in a jump-diffusion model. Dong et al. [Reference Dong, Tang and Wang7] investigated the pricing of vulnerable basket spread options with stochastic liquidity risk. Wang [Reference Wang21] obtained a pricing formula for spread options with stochastically correlated underlying assets.
As is well known, the behavior of real market volatility can not agree with the assumption of constant volatility and thus the above-mentioned formulas have already their own limitation. Refer to, for example, Gatheral [Reference Gatheral11] and Fouque et al. [Reference Fouque, Papanicolaou, Sircar and Sølna9] for empirical evidence of the presence of stochastic factors in the volatility. Furthermore, the prices of many commodity products display seasonality and mean-reversion. So, it is desirable to model the underlying spot prices with stochastic volatility by using another stochastic processes than the geometric Brownian motion even if the bivariate log-normal mixture model of Alexander and Scourse [Reference Alexander and Scourse1] is more elaborate than the earlier standard log-normal scheme. Hikspoors and Jaimungal [Reference Hikspoors and Jaimungal13] and Hikspoors and Jaimungal [Reference Hikspoors and Jaimungal12] used a stochastic volatility model to derive some analytic formula for the commodity spread option prices. Carmona and Sun [Reference Carmona and Sun6] suggested a multiscale stochastic volatility model for the pricing of spread options. In their model, the underlying asset prices are spot prices with linearly growing drift. In most practical cases, however, the underlying asset prices of commodity derivatives are futures prices, not spot prices. Also, it is only futures prices that one can observe in commodity markets.
In this paper, we consider spread options on futures and assume that the underlying asset prices of the futures are mean-reverting and possess a seasonality factor. In addition, the volatilities of the underlying assets are assumed to have both fast and slow variation factors. Approximation approach associated with multiscale stochastic volatility used in this paper is well established in the literature. However, most of the relevant works have been concerned with (vanilla and exotic) options on the spot (index), not on futures. In practice, the mean-reversion (in general, mixing) property, which is the heart of the multiscale stochastic volatility framework, is more important in the price movement of commodity products than securities. So, we are particularly concerned with associating the multiscale stochastic volatility with spread derivatives on futures and come up with this study. We use asymptotic analysis to obtain a closed-form solution for the approximate prices of spread options. As a result, the option price calculation is expected to be much faster than the result of Monte–Carlo simulation and the work of Fouque et al. [Reference Fouque, Saporito and Zubelli10] for single asset options is generalized to multi-asset options and the formula obtained by Kim and Park [Reference Kim and Park14] for exchange options is extended to spread options. These are the main contributions of this paper to the relevant literature. Caldana and Fusai [Reference Caldana and Fusai4] also proposed an exact formula for the approximate spread option prices. Their formula is available only in the case that the joint characteristic function of two underlying log prices is available in closed form like the prices described by the Cox–Ingersoll–Ross (CIR) process. Our formula in this paper is given under a more general condition in the sense that the joint characteristic function is not required and the volatility of the log prices is driven by two different (fast and slow) scales of variations.
The rest of this paper is organized as follows. In Section 2, we formulate the dynamics of underlying asset prices and their stochastic volatilities. In Section 3, the futures prices of these underlying assets are obtained by using an asymptotic expansion method. In addition, the differential forms of the futures prices are obtained as functions of the futures prices (not the spot prices). In Section 4, we extend the formula in Bjerksund and Stensland [Reference Bjerksund and Stensland2] and obtain a closed-form approximation formula under multiscale stochastic volatility. In Section 5, we discuss the Greeks and how to hedge risks under our pricing model and give a brief comment about possible practical applications of our model. In Section 6, the validity of our analytic formula is verified by numerical experiments and our formula is compared with the formula in Bjerksund and Stensland [Reference Bjerksund and Stensland2]. Also, the effects of the scale parameters on the spread option prices are discussed. Finally, Section 6 concludes.
2. Model formulation
The futures contract on an asset is a standardized contract in which both parties agree to trade the asset at future time 
$T$, for a predetermined price, called strike. The future price 
$F_{t,T}$ at time 
$t$ with maturity 
$T$ of the asset is defined as the strike of the futures contract such that no premium is paid at time 
$t$. If 
$S_t$ denotes the asset price at time 
$t$, then the futures price 
$F_{t,T}$ is given by
\begin{equation} F_{t,T}=\mathbb{E}^{Q}[S_T\,|\,\mathcal{F}_t], \end{equation}
where 
$Q$ is a risk-neutral probability measure and 
$\mathcal {F}_t$ denotes a given filtration. Note that the futures price 
$F_{t,T}$ becomes the spot price 
$S_T$ when 
$t=T$, i.e., 
$F_{T,T}=S_T$.
Since we are interested in the pricing of spread options on futures, we use functions 
$F^{(1)}$ and 
$F^{(2)}$ for the prices of two futures. Then the spread (call) option price with option strike 
$K$ and maturity 
$T$ is defined as
\begin{equation} \mathbb{E}^{Q}[ e^{{-}r(T-t)} \max\{F^{(1)}-F^{(2)}-K, 0\}\,|\,\mathcal{F}_t ]. \end{equation}
Note that if 
$F^{(2)}=0$, the spread option collapses into a standard vanilla (call) option on futures whose value has the well-known Black [Reference Black3] formula and if 
$K = 0$, the option becomes into an exchange option whose value is given by the Margrabe [Reference Margrabe18] formula.
We assume that two asset values, 
$S^{(1)}_t$ and 
$S^{(2)}_t$, follow the risk-neutral dynamics given by
\begin{equation} \begin{aligned} S^{(1)}_t & =\exp\{s^{(1)}_t+U^{(1)}_t \},\\ S^{(2)}_t & =\exp\{s^{(2)}_t+U^{(2)}_t \}, \end{aligned} \end{equation}
respectively, where 
$s^{(1)}_t$ and 
$s^{(2)}_t$ are deterministic functions representing seasonality factors and 
$U^{(1)}_t$ and 
$U^{(2)}_t$ are mean-reverting processes whose volatility is driven by two different (fast and slow) scale factors 
$Y_t$ and 
$Z_t$. The following stochastic differential equations (SDEs) represent the precise dynamics of 
$U^{(1)}_t$ and 
$U^{(2)}_t$, respectively.
\begin{equation} \left\{\begin{aligned} dU^{(1)}_t & =\kappa_1(m_1-U^{(1)}_t)\,dt+\eta(Y_t, Z_t)\,dW^{(1)}_t,\\ dU^{(2)}_t & =\kappa_2(m_2-U^{(2)}_t)\,dt+\xi(Y_t, Z_t)\,dW^{(2)}_t,\\ dY_t & =\frac{1}{\epsilon}\alpha(Y_t)\,dt+\frac{1}{\sqrt{\epsilon}}\beta(Y_t)\,dW^{y}_t,\\ dZ_t & =\delta c(Z_t)\,dt+\sqrt{\delta}g(Z_t)\,dW^{z}_t, \end{aligned}\right. \end{equation}
where 
$\kappa _1$, 
$\kappa _2$, 
$m_1$ and 
$m_2$ are constants and 
$\epsilon$ and 
$\delta$ are small positive parameters and 
$W^{(1)}_t$, 
$W^{(2)}_t$, 
$W^{y}_t$ and 
$W^{z}_t$ are standard Brownian motions under the risk-neutral measure 
$\mathbb {Q}$ with correlation structure given by
\begin{equation} \begin{aligned} dW^{(1)}_t\,dW^{(2)}_t & =\rho_{12}\,dt,\quad dW^{(1)}_t\,dW^{y}_t=\rho_{1y}\,dt, \quad dW^{(1)}_t\,dW^{z}_t & =\rho_{1z}\,dt,\\ dW^{(2)}_t\,dW^{y}_t & =\rho_{2y}\,dt, \quad dW^{(2)}_t\,dW^{z}_t=\rho_{2z}\,dt,\quad dW^{y}_t\,dW^{z}_t & =\rho_{yz}\,dt, \end{aligned} \end{equation}
which is assumed to be positive definite. Apart from the existence and uniqueness of the above SDEs, we need the following model assumptions. The functions 
$\alpha (y)$ and 
$\beta (y)$ are assumed to give a guarantee of the existence of a unique invariant distribution 
$\Phi$ of the process 
$Y_t$. The functions 
$c(z)$ and 
$g(z)$ are smooth and at most linearly growing. The functions 
$\eta (y,z)$ and 
$\xi (y, z)$ are positive, bounded and bounded away from zero and they are square-integrable in 
$y$ with respect to the invariant distribution. For simplicity, the market prices of volatility risk are not considered here.
There is a recent study by Schneider and Tavin [Reference Schneider and Tavin20] that also deals with modeling the stochastic volatility of futures prices incorporating a seasonal component. They use a CIR process for the multi-factor variance of the futures prices where the mean-reversion levels of the variance processes represent time-dependent seasonal components. Then the risk-neutral dynamics of the spot price process implied by the futures-based model follow. On the contrary, we use a deterministic function representing a seasonality factor plus a mean-reverting Ornstein–Uhlenbeck(OU) process with multiscale volatility to create a spot pricing model for each of multiple assets and then the futures price is given by the risk-neutral conditional expectation of the terminal spot price.
3. Valuation of futures and risk-neutral dynamics
3.1. Valuation of futures
In this section, we obtain a first-order pricing approximation of the futures of the underlying assets whose values are given by (3).
Let the futures price 
$F^{(i)}$ (
$i=1,2$) of the 
$i$th asset with maturity 
$T$ be denoted by
\begin{equation} F^{(i)}=F^{(i)}(t, U_t^{(1)}, U_t^{(2)}, Y_t, Z_t). \end{equation}
Then from the Markov property we have an expression
\begin{equation} F^{(i)}(t, u_1,u_2,y,z)=\mathbb{E}^{Q}[ S^{(i)}_T \,|\,U_t^{(1)}=u_1, U_t^{(2)}=u_2, Y_t=y, Z_t=z]. \end{equation}
Using the Feyman–Kac theorem, see, for instance, Oksendal [Reference Oksendal19], the futures price 
$F^{(i)}$ satisfies the following partial differential equation (PDE) problem:
\begin{equation} \begin{aligned} & \mathcal{L}^F F^{(i)}(t, u_1, u_2, y, z)=0,\\ & F^{(i)}(T_i, u_1, u_2, y,z)=\exp\{s^{(i)}_{T}+u_i\} \end{aligned} \end{equation}
for 
$i=1, 2$. Here, the operator 
$\mathcal {L}^F$ is given by
\begin{equation} \mathcal{L}^F=\frac{1}{\epsilon}\mathcal{L}^F_0+\frac{1}{\sqrt{\epsilon}}\mathcal{L}^F_1+\mathcal{L}^F_2 +\sqrt{\delta}\mathcal{A}^F_3+\delta\mathcal{A}^F_0+\sqrt{\frac{\delta}{\epsilon}}\mathcal{M}^F_1, \end{equation}
where
\begin{align} \mathcal{L}^F_0 & =\frac{1}{2}\beta^2(y)\frac{\partial^2}{\partial y^2}+\alpha(y)\frac{\partial}{\partial y},\notag\\ \mathcal{L}^F_1 & =\rho_{1y}\beta(y)\eta(y, z)\frac{\partial^2}{\partial u_1 \partial y}+\rho_{2y}\beta(y)\xi(y, z)\frac{\partial^2}{\partial u_2 \partial y},\notag\\ \mathcal{L}^F_2 & =\frac{\partial}{\partial t}+\kappa_1(m_1-u_1)\frac{\partial}{\partial u_1}+\kappa_2(m_2-u_2)\frac{\partial}{\partial u_2},\notag\\ & \quad +\frac{1}{2}\eta^2(y, z)\frac{\partial^2}{\partial u_1^2}+\frac{1}{2}\xi^2(y, z)\frac{\partial^2}{\partial u_2^2}+\rho_{12}\eta(y, z)\xi(y, z)\frac{\partial^2}{\partial u_1 \partial u_2},\\ \mathcal{A}^F_3 & =\rho_{12}g(z)\eta(y, z)\frac{\partial^2}{\partial u_1 \partial z}+\rho_{12}g(z)\xi(y, z)\frac{\partial^2}{\partial u_2 \partial z}, \notag\\ \mathcal{A}^F_0 & =\frac{1}{2}g^2(z)\frac{\partial^2}{\partial z^2}+c(z)\frac{\partial}{\partial z},\notag\\ \mathcal{M}^F_1 & =\rho_{yz}\beta(y)g(z)\frac{\partial^2}{\partial y \partial z}.\notag \end{align} We are interested in an asymptotic expansion of 
$F^{(i)}$ in powers of 
$\epsilon$ and 
$\delta$ in the following form:
\begin{equation} F^{(i)}= \sum_{i, j \geq 0} (\sqrt{\epsilon})^i(\sqrt{\delta})^j F^{(i)}_{ij}. \end{equation}
Then one can expand 
$\mathcal {L}^FF^{(i)}$ as
\begin{align} \mathcal{L}^FF^{(i)}& =\frac{1}{\epsilon}(\mathcal{L}^F_0F_0) +\frac{1}{\sqrt{\epsilon}}(\mathcal{L}^F_0F_{10}+\mathcal{L}^F_1F_0) +(\mathcal{L}^F_0F_{20}+\mathcal{L}^F_1F_{10}+\mathcal{L}^F_2F_{0})\nonumber\\ & \quad +\sqrt{\delta}(\mathcal{L}^F_0F_{21}+\mathcal{L}^F_1F_{11}+\mathcal{L}^F_2F_{01} +\mathcal{A}^F_3F_0+\mathcal{M}^F_1F_{10})\nonumber\\ & \quad +\delta(\mathcal{L}^F_0F_{22}+\mathcal{L}^F_1F_{12}+\mathcal{L}^F_2F_{02}+\mathcal{A}^F_0F_0 +\mathcal{A}^F_3F_{01}+\mathcal{M}^F_1F_{11})\nonumber\\ & \quad +\sqrt{\frac{\delta}{\epsilon}}(\mathcal{L}^F_0F_{11}+\mathcal{L}^F_1F_{01}+\mathcal{M}^F_1F_0) +\sqrt{\epsilon}(\mathcal{L}^F_0F_{30}+\mathcal{L}^F_1F_{20}+\mathcal{L}^F_2F_{10})\nonumber\\ & \quad +\epsilon(\mathcal{L}^F_0F_{40}+\mathcal{L}^F_1F_{30}+\mathcal{L}^F_2F_{20}) +\frac{\sqrt{\delta}}{\epsilon}(\mathcal{L}^F_0F_{01})+\frac{\delta}{\epsilon}(\mathcal{L}^F_0F_{02})\nonumber\\ & \quad +\frac{\delta}{\sqrt{\epsilon}}(\mathcal{L}^F_0F_{12}+\mathcal{L}^F_1F_{02}+\mathcal{M}^F_1F_{01})\nonumber\\ & \quad +\sqrt{\epsilon \delta}(\mathcal{L}^F_0F_{31}+\mathcal{L}^F_1F_{21}+\mathcal{L}^F_2F_{11}+\mathcal{A}^F_3F_{10}+\mathcal{M}^F_1F_{21})+\cdots, \end{align}
where 
$F_{0}$ is used for 
$F_{00}$ and the superscript 
$(i)$ of 
$F_{ij}$ is omitted for simplicity.
If we use the asymptotic analysis of Fouque et al. [Reference Fouque, Papanicolaou, Sircar and Sølna9] with notation 
$\langle \,\cdot \,\rangle$ meaning integration with respect to the invariant distribution 
$\Phi$ of the fast scale process 
$Y_t$, then the leading term 
$F^{(i)}_0$ and the first-order correction terms 
$F^{(i)}_{10}$ and 
$F^{(i)}_{01}$ satisfy the PDE problems
\begin{equation} \begin{aligned} \langle \mathcal{L}^F_2 \rangle F^{(i)}_0(t, u_1, u_2, z) & =0, \quad 0\leq t \lt T, \\ F^{(i)}_0(T,u_1, u_2, z) & =\exp \{ s^{i}_T+u_i \}, \end{aligned} \end{equation}
and
\begin{equation} \begin{aligned} \langle \mathcal{L}^F_2 \rangle \tilde{F}^{(i)}_{10}(t, u_1, u_2, z) & ={-}\langle \mathcal{L}^F_1F^{(i)}_{20}\rangle, \quad 0\leq t \lt T,\\ \tilde{F}^{(i)}_{10}(T, u_1, u_2, z) & =0,\\ \langle \mathcal{L}^F_2 \rangle \tilde{F}^{(i)}_{01}(t, u_1, u_2, z) & ={-}\langle \mathcal{A}^F_3F^{(i)}_{0} \rangle, \quad 0\leq t \lt T,\\ \tilde{F}^{(i)}_{01}(T, u_1, u_2, z) & =0, \end{aligned} \end{equation}
respectively, where 
$\tilde {F}^{(i)}_{10}=\sqrt {\epsilon }F^{(i)}_{10}$ and 
$\tilde {F}^{(i)}_{01}=\sqrt {\delta }F^{(i)}_{01}$. Solutions of these three PDE problems are given by
\begin{equation} \begin{aligned} F^{(1)}_0 & =\exp \left\{s(T)+m_1+(u_1-m_1)e^{-\kappa_1(T-t)}+\frac{\bar{\eta}^2(z)}{4\kappa_1}(1-e^{{-}2\kappa_1(T-t)})\right\}, \\ F^{(2)}_0 & =\exp \left\{s(T)+m_2+(u_2-m_2)e^{-\kappa_2(T-t)}+\frac{\bar{\xi}^2(z)}{4\kappa_2}(1-e^{{-}2\kappa_2(T-t)})\right\}, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} \tilde{F}^{(1)}_{10} & ={-}V^\eta_1(z)\frac{1}{3\kappa_1}(1-e^{{-}3\kappa_1(T-t)})F^{(1)}_0,\\ \tilde{F}^{(2)}_{10} & ={-}V^\xi_2(z)\frac{1}{3\kappa_2}(1-e^{{-}3\kappa_2(T-t)})F^{(2)}_0, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} \tilde{F}^{(1)}_{01} & ={-}\left(\frac{1-e^{-\kappa_1(T-t)}}{2\kappa_1^2} -\frac{1-e^{3\kappa_1(T-t)}}{6\kappa_1^2}\right)\bar{\eta}(z)\bar{\eta}'(z)W^\eta(z)F^{(1)}_0,\\ \tilde{F}^{(2)}_{01} & ={-}\left(\frac{1-e^{-\kappa_2(T-t)}}{2\kappa_2^2} -\frac{1-e^{3\kappa_2(T-t)}}{6\kappa_2^2}\right)\bar{\xi}(z)\bar{\xi}'(z)W^\xi(z)F^{(2)}_0, \end{aligned} \end{equation}
respectively, where 
$\bar {\eta }(z)=\langle \eta (\cdot,z)\rangle$, 
$\bar {\xi }(z)=\langle \xi (\cdot,z)\rangle$ and
\begin{equation} \begin{aligned} V^{\eta}_1(z) & =\tfrac{1}{2}\sqrt{\epsilon}\rho_{1y}\langle\beta\eta\psi^\prime_1\rangle, \quad V^{\xi}_2(z)=\tfrac{1}{2}\sqrt{\epsilon}\rho_{2y}\langle\beta\xi\psi^\prime_2\rangle,\\ W^{\eta}(z) & ={-}\sqrt{\delta}\rho_{1z}g(z)\langle\eta\rangle, \quad W^{\xi}(z)-\sqrt{\delta}\rho_{2z}g(z)\langle\xi\rangle. \end{aligned} \end{equation}
Here, 
$\psi _i$ are the solutions of the Poisson equations
\begin{equation} \begin{aligned} \mathcal{L}_0\psi_1 & =\eta^2-\langle\eta^2\rangle,\\ \mathcal{L}_0\psi_2 & =\xi^2-\langle\xi^2\rangle,\\ \mathcal{L}_0\psi_3 & =\eta\xi-\langle\eta\xi\rangle, \end{aligned} \end{equation}
respectively.
To sum up the above results, we have the following approximation result for the futures prices:
\begin{equation} \begin{aligned} F^{(1)} & \approx F^{(1)}_0+\tilde{F}^{(1)}_{10}+\tilde{F}^{(1)}_{01} =F^{(1)}_0+\sqrt{\epsilon}F^{(1)}_{10}+\sqrt{\delta}F^{(1)}_{01},\\ F^{(2)} & \approx F^{(2)}_0+\tilde{F}^{(2)}_{10}+\tilde{F}^{(2)}_{01} =F^{(2)}_0+\sqrt{\epsilon}F^{(2)}_{10}+\sqrt{\delta}F^{(2)}_{01}. \end{aligned} \end{equation}
3.2. The risk-neutral dynamics of futures prices
Since the futures prices are used for the settlement of spread options, we need to find the risk-neutral dynamics of the futures price 
$F_{t,T}$ to price a given contingent claim. In this section, we present the dynamics of 
$F^{(i)}$ in an SDE form for 
$i=1,2$.
Since 
$F^{(i)}$ is a martingale under the risk-neutral probability measure 
$Q$, its differential should not have a drift term. So, by applying Ito's formula to 
$F^{(i)}$, we have
\begin{equation} \begin{aligned} dF^{(1)} & =\eta(Y_t, Z_t)\frac{\partial F^{(1)}}{\partial u_1}\,dW^{(1)}_t+\frac{1}{\sqrt{\epsilon}}\beta(Y_t)\frac{\partial F^{(1)}}{\partial y}\,dW^y_t+\sqrt{\delta}g(Z_t)\frac{\partial F^{(1)}}{\partial z}\,dW^z_t,\\ dF^{(2)} & =\xi(Y_t,Z_t)\frac{\partial F^{(2)}}{\partial u_2}\,dW^{(2)}_t+\frac{1}{\sqrt{\epsilon}}\beta(Y_t)\frac{\partial F^{(2)}}{\partial y}\,dW^y_t+\sqrt{\delta}g(Z_t)\frac{\partial F^{(2)}}{\partial z}\,dW^z_t. \end{aligned} \end{equation}
In this expression, the partial derivatives are functions of 
$(t, U_t^{1}, U_t^{2}, Y_t, Z_t)$. The goal of this section is to change the spot prices 
$U^{(i)}$ into the futures prices 
$F^{(i)}$.
Let 
$G^{(1)}(t, \cdot, u_2, y, z)$ and 
$G^{(2)}(t, u_1, \cdot, y, z)$ be the inverse functions of 
$F^{(1)}(t, \cdot, u_2, y, z)$ and 
$F^{(2)}(t, u_1, \cdot, y, z)$, respectively. Namely,
\begin{equation} F^{(1)}(t, G^{(1)}(t, \cdot, u_2, y, z), u_2, y, z)=I({\cdot}), \quad F^{(2)}(t, u_1, G^{(2)}(t, u_1, \cdot, y, z), y, z)=I({\cdot}), \end{equation}
where 
$I$ is identity function. As done for 
$F^{(i)}$, we expand 
$G^{(i)}$ also as
\begin{equation} G^{(i)}= \sum_{j, k \geq 0} (\sqrt{\epsilon})^j(\sqrt{\delta})^k G^{(i)}_{jk}, \end{equation}
where 
$G^{(i)}_{00}=G^{(i)}_0$ is chosen to be the inverse function of 
$F^{(i)}_0$ as follows.
\begin{equation} G^{(1)}_0(t, \cdot, u_2, y, z)=( F^{(1)}_0(t, \cdot, u_2, y, z))^{{-}1}, \quad G^{(2)}_0(t, u_1, \cdot, y, z)=( F^{(2)}_0(t, u_1, \cdot, y, z))^{{-}1} \end{equation}
Then the following identities can be easily derived.
\begin{equation} \begin{aligned} G^{(1)}_{10}(t,u_1,u_2,y,z) & ={-}\frac{F^{(1)}_{10}(t, G^{(1)}_0(t, u_1, u_2, y, z), u_2, y, z)}{\frac{\partial F^{(1)}_0}{\partial u_1}(t, G^{(1)}_0(t, u_1, u_2, y, z), u_2, y,z)},\\ G^{(1)}_{01}(t,u_1,u_2,y,z) & ={-}\frac{F^{(1)}_{01}(t, G^{(1)}_0(t, u_1, u_2, y, z), u_2, y, z)}{\frac{\partial F^{(1)}_0}{\partial u_1}(t, G^{(1)}_0(t, u_1, u_2, y, z), u_2, y,z)}. \end{aligned} \end{equation}
\begin{equation} \begin{aligned} G^{(2)}_{10}(t,u_1,u_2,y,z) & ={-}\frac{F^{(2)}_{10}(t, u_1, G^{(2)}_0(t, u_1, u_2, y, z), y, z)}{\frac{\partial F^{(2)}_0}{\partial u_2}(t, u_1, G^{(2)}_0(t, u_1, u_2, y, z),y,z)},\\ G^{(2)}_{01}(t,u_1,u_2,y,z) & ={-}\frac{F^{(2)}_{01}(t, u_1, G^{(2)}_0(t, u_1, u_2, y, z), y, z)}{\frac{\partial F^{(2)}_0}{\partial u_2}(t, u_1, G^{(2)}_0(t, u_1, u_2, y, z), y,z)}. \end{aligned} \end{equation}
Based on the relations
\begin{equation} \begin{aligned} G^{(1)}(t, F^{(1)}(t,u_1,u_2,y,z), u_2, y, z) & =U^{(1)}(t,u_1,u_2,y,z),\\ G^{(2)}(t, u_1,F^{(2)}(t,u_1,u_2,y,z), y, z) & =U^{(2)}(t,u_1,u_2,y,z), \end{aligned} \end{equation}
we now define 
$\phi ^{(i)}_1$, 
$\phi ^{(i)}_2$ and 
$\phi ^{(i)}_3$ as
\begin{equation} \begin{aligned} \phi^{(1)}_1(t, u_1, u_2, y, z) & =\frac{\partial F^{(1)}}{\partial u_1}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z),\\ \phi^{(2)}_1(t, u_1, u_2, y, z) & =\frac{\partial F^{(2)}}{\partial u_2}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z),\\ \phi^{(1)}_2(t, u_1,u_2, y, z) & =\frac{\partial F^{(1)}}{\partial y}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z), \\ \phi^{(2)}_2(t, u_1,u_2, y, z) & =\frac{\partial F^{(2)}}{\partial y}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z), \\ \phi^{(1)}_3(t, u_1,u_2, y, z) & =\frac{\partial F^{(1)}}{\partial z}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z), \\ \phi^{(2)}_3(t, u_1,u_2, y, z) & =\frac{\partial F^{(2)}}{\partial z}(t, G^{(1)}(t, u_1, u_2, y, z), G^{(2)}(t, u_1, u_2, y, z), y, z), \end{aligned} \end{equation}
respectively. Then, we obtain the following desired SDEs for 
$F^{(i)}$:
\begin{equation} \left\{\begin{aligned} dF^{(1)} & =\phi^{(1)}_1(t, F^{(1)}, F^{(2)}, Y_t, Z_t)\eta(Y_t, Z_t)\,dW^1_t & \\ & \quad +\frac{1}{\sqrt{\epsilon}}\phi^{(1)}_2(t, F^{(1)}, F^{(2)}, Y_t, Z_t)\beta(Y_t)\,dW^y_t+\sqrt{\delta}\phi^{(1)}_3(t, F^{(1)}, F^{(2)}, Y_t, Z_t)g(Z_t)\,dW^z_t,\\ dF^{(2)} & =\phi^{(2)}_1(t, F^{(1)}, F^{(2)}, Y_t, Z_t)\xi(Y_t, Z_t)\,dW^2_t\\ & \quad +\frac{1}{\sqrt{\epsilon}}\phi^{(2)}_2(t, F^{(1)}, F^{(2)}, Y_t, Z_t)\beta(Y_t)\,dW^y_t+\sqrt{\delta}\phi^{(2)}_3(t, F^{(1)}, F^{(2)}, Y_t, Z_t)g(Z_t)\,dW^z_t,\\ dY_t & =\frac{1}{\epsilon}\alpha(Y_t)\,dt+\frac{1}{\sqrt{\epsilon}}\beta(Y_t)\,dW^{y}_t,\\ dZ_t & =\delta c(Z_t)\,dt+\sqrt{\delta}g(Z_t)\,dW^{z}_t. \end{aligned}\right. \end{equation}
4. Pricing spread options on futures
In this section, we consider a spread option with maturity 
$T_0$ on two futures contracts whose prices are denoted by 
$F^{(1)}_{t,T}$ and 
$F^{(2)}_{t,T}$, where 
$T_0 \lt T$ and calculate the no-arbitrage price of this option given by
\begin{equation} P(t, x_1, x_2, y, z, T)=e^{{-}r(T_0-t)}\mathbb{E}^{Q}[(F^{(1)}_{T_0,T}-F^{(2)}_{T_0,T}-K)^{+}\,|\, F^{(1)}_{t,T}=x_1, F^{(2)}_{t,T}=x_2, Y_t=y, Z_t=z]. \end{equation}
4.1. PDE for option price
First, the Feynman–Kac representation of this conditional expectation is given by a PDE for 
$P$ satisfying
\begin{align} & \frac{1}{\epsilon}\left[\frac{1}{2}(\phi^{(1)}_2)^2\beta^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}(\phi^{(2)}_2)^2\beta^2\frac{\partial^2}{\partial x_2^2}+\frac{1}{2}\beta^2\frac{\partial^2}{\partial y^2}+\alpha\frac{\partial}{\partial y}\right.\nonumber\\ & \quad \left.+\phi^{(1)}_2\phi^{(2)}_2\beta^2\frac{\partial^2}{\partial x_1 \partial x_2}+\phi^{(1)}_2\beta^2\frac{\partial^2}{\partial x_1 \partial y}+\phi^{(2)}_2\beta^2\frac{\partial^2}{\partial x_2 \partial y}\right]P\nonumber\\ & \quad+\frac{1}{\sqrt{\epsilon}}\left[\rho_{1y}\phi^{(1)}_1\phi^{(1)}_2\eta\beta\frac{\partial^2}{\partial x_1^2}+\rho_{2y}\phi^{(2)}_1\phi^{(2)}_2\xi\beta\frac{\partial^2}{\partial x_2^2}+\rho_{1y}\phi^{(1)}_1\phi^{(2)}_2\eta\beta\frac{\partial^2}{\partial x_1 \partial x_2}\right.\nonumber\\ & \quad \left.+\rho_{2y}\phi^{(1)}_2\phi^{(2)}_1\xi\beta\frac{\partial^2}{\partial x_1 \partial x_2}+\rho_{1y}\phi^{(1)}_1\eta\beta\frac{\partial^2}{\partial x_1 \partial y}+\rho_{2y}\phi^{(2)}_1\xi\beta\frac{\partial^2}{\partial x_2 \partial y}\right]P\nonumber\\ & \quad +\left[\frac{1}{2}(\phi^{(1)}_1)^2\eta^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}(\phi^{(2)}_1)^2\xi^2\frac{\partial^2}{\partial x_2^2}+\rho_{12}\phi^{(1)}_1\phi^{(2)}_1\eta\xi\frac{\partial^2}{\partial x_1 \partial x_2}+\frac{\partial}{\partial t}-r\cdot\right]P\nonumber\\ & \quad +\sqrt{\delta}\left[2\rho_{1z}\phi^{(1)}_1\phi^{(1)}_3\eta g\frac{\partial^2}{\partial x_1^2}+2\rho_{2z}\phi^{(2)}_1\phi^{(2)}_3\xi g\frac{\partial^2}{\partial x_2^2}+\rho_{1z}\phi^{(1)}_1\phi^{(2)}_3\eta g\frac{\partial^2}{\partial x_1 \partial x_2}\right.\nonumber\\ & \quad \left.+\rho_{2z}\phi^{(1)}_3\phi^{(2)}_1\xi g\frac{\partial^2}{\partial x_1 \partial x_2} +\rho_{1z}\phi^{(1)}_1\eta g\frac{\partial^2}{\partial x_1 \partial z}+\rho_{2z}\phi^{(2)}_1\xi g\frac{\partial^2}{\partial x_2 \partial z}\right]P\nonumber\\ & \quad +\delta\left[\frac{1}{2}(\phi^{(1)}_3)^2g^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}(\phi^{(2)}_3)^2g^2\frac{\partial^2}{\partial x_2^2}+\frac{1}{2}g^2\frac{\partial^2}{\partial z^2}+c\frac{\partial}{\partial z}\right.\nonumber\\ & \quad\left. +\phi^{(1)}_3\phi^{(2)}_3g^2\frac{\partial^2}{\partial x_1 \partial x_2}+\phi^{(1)}_3g^2\frac{\partial^2}{\partial x_1 \partial z}+\phi^{(2)}_3g^2\frac{\partial^2}{\partial x_2 \partial z}\right]P \nonumber\\ & \quad +\sqrt{\frac{\delta}{\epsilon}}\left[\rho_{yz}\phi^{(1)}_2\phi^{(1)}_3\beta g\frac{\partial^2}{\partial x_1^2}+\rho_{yz}\phi^{(2)}_2\phi^{(2)}_3\beta g\frac{\partial^2}{\partial x_2^2}+\rho_{yz}\phi^{(1)}_2\phi^{(2)}_3\beta g\frac{\partial^2}{\partial x_1 \partial x_2}\right.\nonumber\\ & \quad +\rho_{yz}\phi^{(1)}_3\phi^{(2)}_2\beta g\frac{\partial^2}{\partial x_1 \partial x_2}+\rho_{yz}\phi^{(1)}_3\beta g\frac{\partial^2}{\partial x_1 \partial y}+\rho_{yz}\phi^{(1)}_2\beta g\frac{\partial^2}{\partial x_1 \partial z}\nonumber\\ & \quad \left.+\rho_{yz}\phi^{(2)}_3\beta g\frac{\partial^2}{\partial x_2 \partial y}+\rho_{yz}\phi^{(2)}_2\beta g\frac{\partial^2}{\partial x_2 \partial z}+\rho_{yz}\beta g\frac{\partial^2}{\partial y \partial z}\right]P=0. \end{align}
Using the expansion of 
$\phi ^{(i)}_j$ (
$i=1,2$ and 
$j=1,2,3$) in order of 
$\sqrt {\epsilon }$ and 
$\sqrt {\delta }$ and discarding 
$\mathcal {O}(\epsilon +\delta )$ terms, we obtain
\begin{equation} \left\{\begin{aligned} & \left(\frac{1}{\epsilon}\mathcal{L}_0+\frac{1}{\sqrt{\epsilon}}\mathcal{L}_1+\mathcal{L}_2 +\sqrt{\epsilon}\mathcal{L}_3+\sqrt{\delta}\mathcal{A}_3+\sqrt{\frac{\delta}{\epsilon}}\mathcal{M}_1 +\frac{\delta}{\sqrt{\epsilon}}\mathcal{M}_2\right)P=0, \quad t \lt T_0,\\ & P(T_0, x_1, x_2, y,z,T)=(x_1-x_2-K)^{+}, \end{aligned}\right. \end{equation}
where
\begin{align} \mathcal{L}_0& =\mathcal{L}^F_0=\frac{1}{2}\beta^2\frac{\partial^2}{\partial y^2}+\alpha\frac{\partial}{\partial y},\nonumber\\ \mathcal{L}_1& =\rho_{1y}\phi^{(1)}_{100}\eta\beta\frac{\partial^2}{\partial x \partial y}+\rho_{2y}\phi^{(2)}_{100}\xi\beta\frac{\partial^2}{\partial x \partial y},\nonumber\\ \mathcal{L}_2& =\frac{1}{2}(\phi^{(1)}_{100})^2\eta^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}(\phi^{(2)}_{100})^2\xi^2\frac{\partial^2}{\partial x_2^2}+\rho_{12}\phi^{(1)}_{100}\phi^{(2)}_{100}\eta\xi\frac{\partial^2}{\partial x_1 \partial x_2}\nonumber\\ & \quad +\phi^{(1)}_{220}\beta^2\frac{\partial^2}{\partial x_1 \partial y}+\phi^{(2)}_{220}\beta^2\frac{\partial^2}{\partial x_2 \partial y}+\frac{\partial}{\partial t}-r\cdot,\nonumber\\ \mathcal{L}_3& =\rho_{1y}\phi^{(1)}_{100}\phi^{(1)}_{220}\eta\beta\frac{\partial^2}{\partial x_1^2}+\rho_{2y}\phi^{(2)}_{100}\phi^{(2)}_{220}\xi\beta\frac{\partial^2}{\partial x_2^2}+(\rho_{1y}\phi^{(1)}_{100}\phi^{(2)}_{220}\eta\beta+\rho_{2y}\phi^{(1)}_{220}\phi^{(2)}_{100}\xi\beta) \frac{\partial^2}{\partial x_1 \partial x_2}\nonumber\\ & \quad +(\phi^{(1)}_{230}\beta^2+\rho_{1y}\phi^{(1)}_{120}\eta\beta)\frac{\partial^2}{\partial x_1 \partial y}+(\phi^{(2)}_{230}\beta^2+\rho_{2y}\phi^{(2)}_{120}\xi\beta)\frac{\partial^2}{\partial x_2 \partial y},\nonumber\\ \mathcal{A}_3& =(\phi^{(1)}_{221}\beta^2+\rho_{1y}\phi^{(1)}_{111}\eta\beta+\rho_{yz}\phi^{(1)}_{310}\beta g)\frac{\partial^2}{\partial x_1 \partial y}+(\phi^{(2)}_{221}\beta^2+\rho_{2y}\phi^{(2)}_{111}\xi\beta+\rho_{yz}\phi^{(2)}_{310}\beta g)\frac{\partial^2}{\partial x_2 \partial y}\nonumber\\ & \quad +\rho_{1z}\phi^{(1)}_{100}\eta g\frac{\partial^2}{\partial x_1 \partial z}+\rho_{2z}\phi^{(2)}_{100}\xi g\frac{\partial^2}{\partial x_2 \partial z}+2\rho_{1z}\phi^{(1)}_{100}\phi^{(1)}_{300}\eta g\frac{\partial^2}{\partial x_1^2}+2\rho_{2z}\phi^{(2)}_{100}\phi^{(2)}_{300}\xi g\frac{\partial^2}{\partial x_2^2}\nonumber\\ & \quad +(\rho_{1z}\phi^{(1)}_{100}\phi^{(2)}_{300}\eta g+\rho_{2z}\phi^{(1)}_{300}\phi^{(2)}_{100}\xi g)\frac{\partial^2}{\partial x_1 \partial x_2}, \nonumber\\ \mathcal{M}_1& =\rho_{yz}\phi^{(1)}_{300}\beta g\frac{\partial^2}{\partial x_1 \partial y}+\rho_{yz}\phi^{(2)}_{300}\beta g\frac{\partial^2}{\partial x_2 \partial y}+\rho_{yz}\beta g\frac{\partial^2}{\partial y \partial z},\nonumber\\ \mathcal{M}_2& =\rho_{yz}\phi^{(1)}_{301}\beta g\frac{\partial^2}{\partial x_1 \partial y}+\rho_{yz}\phi^{(2)}_{301}\beta g\frac{\partial^2}{\partial x_2 \partial y}. \end{align}
Here, the expansions 
$\phi ^{(i)}_1=\phi ^{(i)}_{100}+\sqrt {\epsilon }\phi ^{(i)}_{110}+\sqrt {\delta }\phi ^{(i)}_{101}+\cdots$, 
$\phi ^{(i)}_2=\phi ^{(i)}_{200}+\sqrt {\epsilon }\phi ^{(i)}_{210}+\sqrt {\delta }\phi ^{(i)}_{201}+\cdots$ and 
$\phi ^{(i)}_3=\phi ^{(i)}_{300}+\sqrt {\epsilon }\phi ^{(i)}_{310}+\sqrt {\delta }\phi ^{(i)}_{301}+\cdots$ have been used. We shall see that some terms in these expansions are not necessary to calculate for a first-order approximation of our interest (because of independence on the variable 
$y$). So, we express explicitly only terms that are going to be used in the following section. They are
\begin{equation} \begin{aligned} \phi^{(1)}_{100} & =e^{-\kappa_1(T-t)}x_1, \quad \phi^{(2)}_{100}=e^{-\kappa_2(T-t)}x_2,\\ \phi^{(1)}_{220} & ={-}\frac{1}{2}e^{{-}2\kappa_1(T-t)}\frac{\partial \psi_1}{\partial y}(y, z)x_1, \quad \phi^{(2)}_{220}={-}\frac{1}{2}e^{{-}2\kappa_2(T-t)}\frac{\partial \psi_2}{\partial y}(y, z)x_2,\\ \phi^{(1)}_{300} & =\frac{1}{2\kappa_1}\bar{\eta}(z)\bar{\eta}'(z)(1-e^{{-}2\kappa_1(T-t)})x_1,\quad \phi^{(2)}_{300}=\frac{1}{2\kappa_2}\bar{\xi}(z)\bar{\xi}'(z)(1-e^{{-}2\kappa_2(T-t)})x_2. \end{aligned} \end{equation}
4.2. First-order approximation
For the spread option price 
$P$, we are interested in an asymptotic expansion given by
\begin{equation} P(t,x_1,x_2,y,z,T)=\sum_{i, j \geq 0} (\sqrt{\epsilon})^i(\sqrt{\delta})^j P_{ij}(t,x_1,x_2,y,z,T). \end{equation}
Substituting this expansion into the PDE (32), we have
\begin{align} & \frac{1}{\epsilon}(\mathcal{L}_0P_{00}) +\frac{1}{\sqrt{\epsilon}}(\mathcal{L}_0P_{10}+\mathcal{L}_1P_{00}) +(\mathcal{L}_0P_{20}+\mathcal{L}_1P_{10}+\mathcal{L}_2P_{00})\nonumber\\ & \quad +\sqrt{\epsilon}(\mathcal{L}_0P_{30}+\mathcal{L}_1P_{20}+\mathcal{L}_2P_{10}+\mathcal{L}_3P_{00})\nonumber\\ & \quad +\sqrt{\delta}(\mathcal{L}_0P_{21}+\mathcal{L}_1P_{11}+\mathcal{L}_2P_{01}+\mathcal{A}_3P_{00}+\mathcal{M}_1P_{10})\nonumber\\ & \quad +\frac{\sqrt{\delta}}{\epsilon}(\mathcal{L}_0P_{01}) +\sqrt{\frac{\delta}{\epsilon}}(\mathcal{L}_0P_{11}+\mathcal{L}_1P_{01}+\mathcal{M}_1P_{00})\nonumber\\ & \quad +\frac{\delta}{\sqrt{\epsilon}}(\mathcal{L}_0P_{12}+\mathcal{L}_1P_{02}+\mathcal{M}_1P_{01} +\mathcal{M}_2P_{00})+\cdots=0. \end{align}
From (36), one can obtain a system of PDEs, 
$\mathcal {L}_0P_{00}=0$, 
$\mathcal {L}_0P_{10}+\mathcal {L}_1P_{00}=0$, 
$\mathcal {L}_0P_{20}+\mathcal {L}_1P_{10}+\mathcal {L}_2P_{00}=0$ and so on. Using the asymptotic analysis of Fouque et al. [Reference Fouque, Papanicolaou, Sircar and Sølna9] for this PDE system, one can find solutions for 
$P_{00}$, 
$P_{10}$ and 
$P_{01}$ as follows.
First, if we assume that 
$P_{00}$ does not grow exponentially in 
$y$, then one can find that 
$P_{00}$ is independent of 
$y$ and 
$P_{00}$ satisfies the PDE problem
\begin{equation} \left\{\begin{aligned} & \langle \mathcal{L}_2 \rangle P_{00}(t,x_1,x_2,z,T)=0, \quad t \lt T_0, \\ & P_{00}(T_0, x_1, x_2,z,T)=\max\{x_1-x_2-K, 0\}, \end{aligned}\right. \end{equation}
where
\begin{align} \langle \mathcal{L}_2 \rangle& =\frac{1}{2}e^{{-}2\kappa_1(T-t)}\bar{\eta}^2x_1^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}e^{{-}2\kappa_2(T-t)}\bar{\xi}^2x_2^2\frac{\partial^2}{\partial x_2^2}\nonumber\\ & \quad +\rho_{12}e^{-(\kappa_1+\kappa_2)(T-t)}\langle \eta\xi \rangle x_1x_2\frac{\partial^2}{\partial x_1 \partial x_2}+\frac{\partial}{\partial t}-r. \end{align}
Then 
$P_{00}$ can be approximated by the spread option price 
$P_0$ of Bjerksund and Stensland [Reference Bjerksund and Stensland2] extending Kirk's formula as shown in Appendix. Using more intuitive notation 
$P_{\mathrm {KBS}}$, it is expressed as
\begin{align} P_0& =P_{\mathrm{KBS}}(t, x_1, x_2, T; \bar{\sigma}_{1,\tau}, \bar{\sigma}_{2,\tau})\nonumber\\ & =e^{{-}r(T_0-t)}(x_1N(d_1)-x_2N(d_2)-KN(d_3)), \end{align}
where 
$N(\cdot )$ denotes the standard normal cumulative probability function and 
$d_1$, 
$d_2$ and 
$d_3$ are given by
\begin{equation} \begin{aligned} & d_1=\frac{\log(x_1/a)+( \frac{1}{2}\bar{\sigma}_{1,\tau}^2-b\rho_{12}\bar{\sigma}_{1,\tau}\bar{\sigma}_{2,\tau}+\frac{1}{2}b^2\bar{\sigma}_{2,\tau}^2)(T_0-t)}{\sigma\sqrt{T_0-t}},\\ & d_2=\frac{\log(x_1/a)+( -\frac{1}{2}\bar{\sigma}_{1,\tau}^2+\rho_{12}\bar{\sigma}_{1,\tau}\bar{\sigma}_{2,\tau}+\frac{1}{2}b^2\bar{\sigma}_{2,\tau}^2-b\bar{\sigma}_{2,\tau}^2)(T_0-t)}{\sigma\sqrt{T_0-t}},\\ & d_3=\frac{\log(x_1/a)+( -\frac{1}{2}\bar{\sigma}_{1,\tau}^2+\frac{1}{2}b^2\bar{\sigma}_{2,\tau}^2)(T_0-t)}{\sigma\sqrt{T_0-t}}, \end{aligned} \end{equation}
respectively, where
\begin{equation} \begin{aligned} & \sigma^2=\bar{\sigma}_{1,\tau}^2+2b\rho_{12}\bar{\sigma}_{1,\tau}\bar{\sigma}_{2,\tau}+b^2\bar{\sigma}_{2,\tau}^2,\\ & a=x_2+K, \\ & b=\frac{x_2}{x_2+K}. \end{aligned} \end{equation}
Here, 
$\bar {\sigma }_{1,\tau }$ and 
$\bar {\sigma }_{2,\tau }$ are defined in Appendix.
Next, by assuming the same growth condition as for 
$P_{00}$ with respect to variable 
$y$, 
$P_{10}$ is also independent of 
$y$. From the above result for 
$P_{00}$ and the PDE 
$\mathcal {L}_0P_{30}+\mathcal {L}_1P_{20}+\mathcal {L}_2P_{10}+\mathcal {L}_3P_0=0$ with the terminal condition 
$P_{10}(T_0, x_1, x_2, z,T)=0$, one can find 
$\tilde {P}_{10}:=\sqrt {\epsilon }P_{10}$ as
\begin{align} \tilde{P}_{10}& =\nu_{3, 0}\mathcal{R}^{\epsilon}_1(z)D^{(1)}_1D^{(1)}_2P_0+\nu_{1, 2}\mathcal{R}^{\epsilon}_2(z)D^{(1)}_1D^{(2)}_2P_0\nonumber\\ & \quad +\nu_{2,1}\mathcal{R}^{\epsilon}_3(z)D^{(1)}_1D^{(1)}_1D^{(2)}_1P_0+\nu_{2, 1}\mathcal{R}^{\epsilon}_4(z)D^{(1)}_2D^{(2)}_1P_0\nonumber\\ & \quad +\nu_{0,3}\mathcal{R}^{\epsilon}_5(z)D^{(2)}_1D^{(2)}_2P_0+\nu_{1, 2}\mathcal{R}^{\epsilon}_6(z)D^{(1)}_1D^{(2)}_1D^{(2)}_1P_0\nonumber\\ & \quad +\nu_{3,0}\mathcal{R}^{\epsilon}_1(z)D^{(1)}_2P_0+\nu_{0, 3}\mathcal{R}^{\epsilon}_5(z)D^{(2)}_2P_0\nonumber\\ & \quad +\nu_{1, 2}\mathcal{R}^{\epsilon}_2(z)D^{(1)}_1D^{(2)}_1P_0+\nu_{2, 1}\mathcal{R}^{\epsilon}_4(z)D^{(1)}_1D^{(2)}_1P_0, \end{align}
where
\begin{equation} \begin{aligned} \nu_{i, j} & =\nu(i\kappa_1+j\kappa_2; t, T_0, T),\\ \nu(\kappa; t, T_0, T) & :=\frac{e^{-\kappa(T-T_0)}-e^{-\kappa(T-t)}}{\kappa}, \end{aligned} \end{equation}
as shown in Appendix. Here, 
$\mathcal {R}^{\epsilon }_{i}$ 
$(i=1,2,3,4,5)$ and 
$D^{(j)}_n$ 
$(j=1,2, n=1,2)$ are given in Appendix, respectively.
By assuming the same growth condition as for 
$P_{00}$ and 
$P_{10}$ with respect to variable 
$y$, 
$P_{01}$ is also independent of 
$y$. The above results for 
$P_0$ and 
$P_{10}$ and the PDE 
$\mathcal {L}_0P_{21}+\mathcal {L}_1P_{11}+\mathcal {L}_2P_{01}+\mathcal {A}_3P_0+\mathcal {M}_1P_{10}=0$ with the terminal condition 
$P_{01}(T_0, x_1, x_2, z,T)=0$ provide us with the solution 
$\tilde {P}_{01}:=\sqrt {\delta }P_{01}$ given by
\begin{align} \tilde{P}_{01}& =\mu_1\mathcal{R}_1^\delta(z)D^{(1)}_1\frac{\partial P_0}{\partial z}+\mu_2\mathcal{R}_2^\delta(z)D^{(2)}_1\frac{\partial P_0}{\partial z}\nonumber\\ & \quad +(\nu_{1,0}-\nu_{3,0})\mathcal{R}_3^\delta(z)D^{(1)}_2P_0+(\nu_{0,1}-\nu_{0,3})\mathcal{R}_4^\delta(z)D^{(2)}_2P_0\nonumber\\ & \quad +(\nu_{1,0}-\nu_{1,2})\mathcal{R}_5^\delta(z)D^{(1)}_1D^{(2)}_1P_0+(\nu_{0,1}-\nu_{2,1})\mathcal{R}_6^\delta(z)D^{(1)}_1D^{(2)}_1P_0, \end{align}
where
\begin{equation} \begin{aligned} \mu_1 & =\frac{(T_0-t)f_{1,0}-\nu_{1,0}}{\kappa_1(T_0-t)},\\ \mu_2 & =\frac{(T_0-t)f_{0,1}-\nu_{0,1}}{\kappa_2(T_0-t)}. \end{aligned} \end{equation}
as shown in Appendix. Here, 
$\mathcal {R}^{\delta }_{i}$ 
$(i=1,2,3,4,5,6)$ and 
$D^{(j)}_n$ 
$(j=1,2, n=1,2)$ are given in Appendix, respectively.
To sum up the above results, we obtain a closed-form approximation for the spread option price 
$P$ expressed by
\begin{equation} P \approx P_0+\tilde{P}_{10}+\tilde{P}_{01}=P_{0}+\sqrt{\epsilon}P_{10}+\sqrt{\delta}P_{01}, \end{equation}
where the leading order term 
$P_{0}$ and the first-order correction terms 
$\tilde {P}_{10}$ and 
$\tilde {P}_{01}$ are given by the formulas (39), (42) and (44), respectively. The approximation 
$P_0+\tilde {P}_{10}+\tilde {P}_{01}$ can be easily calculated by taking derivatives of 
$P_{0}$ (the KBS formula) with respect to 
$x_1$, 
$x_2$ and 
$z$.
Remark. Since the pricing formula (46) is not an exact solution but it is a first-order approximation result derived formally, an accuracy analysis for the error is desirable. However, a rigorous argument would follow the line of Fouque et al. [Reference Fouque, Papanicolaou, Sircar and Solna8] or Appendix B of Fouque et al. [Reference Fouque, Saporito and Zubelli10] except that two underlying assets require a little bit more calculation than the one asset case. Also, the accuracy will be confirmed indirectly through Monte–Carlo simulation as shown in the next section. So, we just state the accuracy result here without proof.
\begin{equation} |P-(P_0+\tilde{P}_{10}+\tilde{P}_{01})|= \mathcal{O}(\epsilon+\delta+\sqrt{\epsilon\delta}). \end{equation}
5. Computation of the Greeks and hedging
For practical use of the pricing formula, we need to calculate the Greeks, i.e., the sensitivity of the prices with respect to underlyings, time and volatilities. For the KBS formula, it is simple to derive closed-form Greeks (see [Reference Bjerksund and Stensland2]). On the other hand, our formula (46) has the form of a multiscale extension with the perturbation terms 
$P_{10}$ and 
$P_{01}$. These perturbation terms are too complicated to calculate their derivatives in closed-form. Therefore, we need to use a finite difference scheme to compute the Greeks. For example, the delta with respect to underlying asset 1 is calculated by the central difference
\begin{equation} \frac{\partial P}{\partial F^{(1)}} \approx \frac{P(F^{(1)}+dF^{(1)})-P(F^{(1)}-dF^{(1)})}{2dF^{(1)}}, \end{equation}
where 
$dF^{(1)}=0.01\times F^{(1)}$. We note that by convention a 1% change in 
$F^{(1)}$ is used by practitioners and in fact gamma cash is defined as the change in delta cash for a 1% move in the underlying. Other Greeks can be computed similarly. In terms of time-efficiency, it is much faster than Monte-Carlo (MC) simulation although it takes more time than the case knowing closed-form Greeks.
Briefly speaking of how to use the sensitivities in hedging procedure, one can think of a hedging strategy of oil companies (or power-plant companies) as follows. These companies have a risk exposure to delta spread since the crude oil and the refined oil prices move differently in the real market. So, these companies could buy a put spread option to hedge the downside risk of their position.
6. Numerical results
In Section 4, we have obtained a first-order approximation formula for the price of a spread option. In this section, we check the validity of the formula and investigate the sensitivity of the option price to some parameters through some numerical experiments and also discuss a calibration procedure based on crude oil market data to estimate the pricing parameters.
6.1. Validity of the formula and price sensitivity
We investigate the sensitivity of the spread option price with respect to the correlation between two underlyings and the mean-reversion speed the fast scale process and verify the accuracy of our analytic formula via MC simulation.
Since the model (4) contains general functions, we need to specify those functions to calculate the formula and utilize the MC simulation. The functions in the model dynamics are chosen as 
$\zeta (y)=\exp {(y)}$, 
$\zeta _1(z)=\zeta _2(z)=\exp {(z)}$, 
$\alpha (y)=m_y-y$, 
$\beta (y)=\nu _y\sqrt {2}$, 
$c(z)=0$, 
$g(z)=z$ based on the specification 
$\eta (y, z)=\zeta (y)\zeta _1(z)$ and 
$\xi (y, z)=\zeta (y)\zeta _2(z)$ in the model and a fast mean-reverting OU diffusion for the fast scale process 
$Y_t$ and a geometric Brownian motion for the slow scale process 
$Z_t$. The model parameters are chosen as 
$m_y=-0.7$, 
$\nu _y=0.5$ and 
$s_1(t)=s_2(t)=0$. In addition, 
$\kappa _1=\kappa _2=1$, 
$m_1=4.5$, 
$m_2=4.4$, 
$\rho _{1y}=\rho _{1z}=\rho _{2y}=\rho _{2z}=-0.5$, 
$\rho _{yz}=0$, 
$U^{(1)}_0=4.8$, 
$U^{(2)}_0=4.7$, 
$Y_0=-0.8$, 
$Z_0=-0.78$. The risk-neutral interest rate is 
$r=0.05$ and 
$T=T_0+30/365$ is chosen since the futures expiry is usually one month after the option expiry. The choice of zero-valued seasonality functions 
$s_1(t)$ and 
$s_2(t)$ is based on the historical WTI and Brent futures data. As shown in Figure 1, there is no seasonality in the futures term-structure for the period from May 2022 to May 2024. One can also see the zeroing of seasonality in Fouque et al. [Reference Fouque, Saporito and Zubelli10].
Figure 1. The term-structures of WTI and Brent crude oil futures for the period from May 2022 to May 2024
For more precise simulation, we use the control variates method, one of a well-known variance reduction technique used in MC simulation. Let the spread option price with constant volatilities be a control variate, since its analytic solution is already known as the KBS formula 
$P_{{\rm KBS}}$ and its payoff is highly correlated with our target payoff driven by stochastic volatility processes. Specifically, let 
$h_{{\rm KBS}}(T_0)$ and 
$h(T_0)$ be the payoffs of the spread options with constant and multiscale stochastic volatility, respectively. Written in an equation, we have
\begin{equation} \begin{aligned} P_{{\rm KBS}} & =e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h_{{\rm KBS}}(T_0)],\\ P & =e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h(T_0)]. \end{aligned} \end{equation}
Then it follows that
\begin{align} P& =e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h(T_0)-h_{{\rm KBS}}(T_0)+h_{{\rm KBS}}(T_0)]\nonumber\\ & =e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h_{{\rm KBS}}(T_0)] + e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h(T_0)-h_{{\rm KBS}}(T_0)]\nonumber\\ & =P_{{\rm KBS}} + e^{{-}r(T_0-t)}\mathbb{E}^{Q}[h(T_0)-h_{{\rm KBS}}(T_0)]. \end{align}
So, we simulate the last expectation term.
Figure 2 shows the spread option prices given by the approximated formula with respect to the correlation 
$\rho _{12}$ between two underlyings and the strike 
$K$ of the spread option. As desired, the spread option price decreases as two underlyings get correlated more positively regardless of the moneyness.
Figure 2. The surface of the spread option price as a function of the correlation between two underlyings and the moneyness with 
$\epsilon =0.1$, 
$\delta =0.01$ and 
$T_0=0.5$.
On the other hand, Table 1 presents the spread option prices from the KBS result, our pricing formula 
$P$ and MC simulation for a number of correlation values between two underlyings. The prices from our formula tend to match well with the MC simulation results compared to the KBS prices, except the highly correlated (
$\rho _{12} \geq 0.95$) or deep-OTM (moneyness 
$\geq 1.4$) cases as shown in Table 1.
Table 1. The spread option prices from the Kirk–Bjerksund–Stensland (KBS) result, our pricing formula and MC simulation with 
$\epsilon =0.1$, 
$\delta =0.01$ and 
$T_0=0.5$.
The spread option price varies with the two scale parameters 
$\epsilon$ and 
$\delta$ of our multiscale volatility model. So, the pricing formula has a higher degree of freedom than the KBS formula. As shown in Figure 3, the option price gets larger regardless of moneyness as the fast scale parameter 
$\epsilon$ (the reciprocal of the speed of mean-reversion) decreases. This is naturally justified by the fact that the volatility is increasing and approaching the mean level which is higher than the initial value as the volatility is more highly mean-reverting. This phenomenon is more noticeable when the two scale factors become more negatively correlated.
Figure 3. The behavior of the spread option prices for various values of 
$\epsilon$ with 
$\delta =0.01$ and 
$T_0=0.5$.
6.2. Calibration to real market data
In this section, we discuss an estimation procedure of our model on real data. The WTI-Brent spread option is not liquid enough to calibrate the model from its historical data. One possible alternative to use is the WTI and Brent crude oil data to estimate the parameters. The WTI and Brent crude oil implied volatilities, shown in the Bloomberg terminal calculated from the option market prices, are chosen for the required data. These volatilities are converted from the historical option prices by using the Black model (not the Black–Scholes model) since the underlyings are futures (not spots). By using vanilla option data, one can get more market information than using the spread option data only. The general validity of applying the vanilla option parameters to exotic options is explained in Fouque et al. [Reference Fouque, Papanicolaou, Sircar and Sølna9] and the detailed procedure can be found in the work of Carmona and Sun [Reference Carmona and Sun6] and Fouque et al. [Reference Fouque, Saporito and Zubelli10].
Instead of estimating each of all model parameters directly, a desirable calibration procedure should aim to obtain the group market parameters such as 
$\mathcal {R}^{\epsilon }_i(z)$'s and 
$\mathcal {R}^{\delta }_i(z)$'s that appear in the correction terms 
$\tilde {P}_{10}$ and 
$\tilde {P}_{01}$. First, we utilize WTI-option volatility data to obtain the parameters 
$\kappa _1$, 
$\bar {\eta }(z)$, 
$\mathcal {R}^{\epsilon }_1(z)$ and 
$\mathcal {R}^{\delta }_3(z)$ of the first asset 
$F^{(1)}$. Then by using Brent-option volatility data, the parameters 
$\kappa _2$, 
$\bar {\xi }(z)$, 
$\mathcal {R}^{\epsilon }_5(z)$ and 
$\mathcal {R}^{\delta }_4(z)$ of the second asset 
$F^{(2)}$ are estimated. The remaining group market parameters can be calculated from their definitions given in Appendix.
To estimate the desired parameters from each vanilla option data, we follow the lines of Fouque et al. [Reference Fouque, Saporito and Zubelli10], in which the calibration procedure consists of two steps. The results of the first step are drawn in Figure 4. The dotted plots represent market implied volatilities with respect to the log-moneyness-to-maturity ratio (LMMR) on May 6, 2022. For each maturity, volatility skews are fitted to the market data. As the second step, we use the parameters obtained at the first step to obtain the final parameters as follows. For WTI-options, 
$\kappa _1=0.830005$, 
$\bar {\eta }(z)=0.667161$, 
$\mathcal {R}^{\epsilon }_1(z)=0.015426$ and 
$\mathcal {R}^{\delta }_3(z)=-0.166520$. For Brent options, 
$\kappa _2=1.381936$, 
$\bar {\xi }(z)=0.761285$, 
$\mathcal {R}^{\epsilon }_5(z)=0.022263$ and 
$\mathcal {R}^{\delta }_4(z)=-0.140858$. By using these group parameters and our analytic formula, we can calculate the spread option prices. The resultant option price surface is shown in Figure 5. Note that the option price increases with the moneyness since the ATM underlying spread value is negative.
Figure 4. WTI and Brent option implied volatilities (dotted lines) and calibrated volatility skews (solid lines) for several choices of maturity.
Figure 5. Spread option prices against moneyness and time to maturity based on the calibrated parameters.
Based on the calibrated parameters, one can test the accuracy of the approximated prices. We obtain the theoretical prices and the real market prices with their errors in Figure 6. As mentioned earlier, we have only a few available price data of the WTI-Brent spread options. All the strike levels are converted to the moneyness. In July 2022 corresponding to a very short time to maturity, the deep out-of-the money options show a big error. However, one can find that the error tends to become seriously smaller as time to maturity gets longer.
Figure 6. Market prices and spread option prices with the calibrated parameters.
7. Conclusion
In this work, we consider the problem of pricing spread options on futures contracts. We assume mean-reverting dynamics for the futures prices with multiscale stochastic volatility. This type of model framework is inspired by the commodity markets, in which the futures contract is a practical choice for the underlying assets of spread options. The analytic pricing result obtained in this article is an improvement of not only the classical Kirk's formula but also its improvement by Bjerksund and Stensland [Reference Bjerksund and Stensland2] in two aspects. First, our approximation formula is still given in a closed-form solution and yet better than their results in terms of accuracy. Second, our pricing formula is about the pricing of a derivative (spread option) of a derivative (futures), which is thought of as a composite contingent claim, instead of the valuation of a derivative of an underlying asset traded in a spot market. Our result is mathematically more accurate than the known results and financially more well suited for application to the commodity markets.
Acknowledgments
The research of J.-H. Kim was supported by the National Research Foundation of Korea NRF2021R1A2C1004080.
Competing interests
The authors declare no conflict of interest.
Appendix
A.1. Derivation of solution 
$P_{00}$
If 
$\langle \mathcal {L}_2 \rangle$ has the form of the operator 
$\mathcal {L}_{\mathrm {bln}}$ corresponding to a bivariate log-normal model, where
\begin{equation} \mathcal{L}_{\mathrm{bln}}(\sigma_1, \sigma_2)=\frac{1}{2}\sigma_1^2x_1^2\frac{\partial^2}{\partial x_1^2}+\frac{1}{2}\sigma_2^2x_2^2\frac{\partial^2}{\partial x_2^2} +\rho\sigma_1\sigma_2x_1x_2\frac{\partial^2}{\partial x_1 \partial x_2}+\frac{\partial}{\partial t}-r, \end{equation}
then it would be easier to obtain a solution for 
$P_{00}$ from the result of Kirk [Reference Kirk15] or Bjerksund and Stensland [Reference Bjerksund and Stensland2]. However, it is not the case since 
$\langle \eta \xi \rangle \neq \bar {\eta }\bar {\xi }$ holds. To handle this problem, we assume that 
$\eta$ and 
$\xi$ are given by the separable form
\begin{equation} \begin{aligned} \eta(y, z) & =\zeta(y)\zeta_1(z),\\ \xi(y, z) & =\zeta(y)\zeta_2(z) \end{aligned} \end{equation}
and we define
\begin{equation} \begin{aligned} \bar{\sigma}_1(t, z) & =e^{-\kappa_1(T-t)}\bar{\zeta}\zeta_1(z),\\ \bar{\sigma}_2(t, z) & =e^{-\kappa_2(T-t)}\bar{\zeta}\zeta_2(z),\\ (\bar{\zeta})^2 & :=\langle \zeta^2 \rangle. \end{aligned} \end{equation}
Note that the volatilities 
$\bar {\sigma }_1(t, z)$ and 
$\bar {\sigma }_2(t, z)$ still have dependence on time 
$t$. In order to apply the formula of Bjerksund and Stensland [Reference Bjerksund and Stensland2], the volatilities should be constant. So, we define the time-averaged volatilities
\begin{equation} \begin{aligned} \bar{\sigma}_{1,\tau}^2 & =\frac{1}{T_0-t}\int_t^{T_0}\bar{\sigma}_1^2(s, z)\,ds,\\ \bar{\sigma}_{2,\tau}^2 & =\frac{1}{T_0-t}\int_t^{T_0}\bar{\sigma}_2^2(s, z)\,ds. \end{aligned} \end{equation}
Then, for each fixed 
$t$, 
$P_{00}$ is approximated by the solution 
$P_0$ of the PDE problem
\begin{equation} \left\{\begin{aligned} & \mathcal{L}_{\mathrm{bln}}(\bar{\sigma}_{1,\tau}, \bar{\sigma}_{2,\tau})P_0(t,x_1,x_2,z,T)=0, \\ & P_0(T_0, x_1, x_2, z,T)=\max\{x_1-x_2-K, 0\} \end{aligned}\right. \end{equation}
whose solution is exactly the same as the spread option price of Bjerksund and Stensland [Reference Bjerksund and Stensland2] extending Kirk's formula and it is expressed as (39).
A.2. Derivation of solution $P_{10}$
From the result for 
$P_{00}$ and the PDE 
$\mathcal {L}_0P_{30}+\mathcal {L}_1P_{20}+\mathcal {L}_2P_{10}+\mathcal {L}_3P_0=0$, one can find that 
$\tilde {P}_{10}:=\sqrt {\epsilon }P_{10}$ satisfies the PDE
\begin{align} \langle \mathcal{L}_2 \rangle \tilde{P}_{10}& ={-}f_{3, 0}(t, T)\mathcal{R}^{\epsilon}_{1}D^{(1)}_1D^{(1)}_2P_0-f_{1, 2}(t, T)\mathcal{R}^{\epsilon}_{2}D^{(1)}_1D^{(2)}_2P_0\nonumber\\ & \quad -f_{2, 1}(t, T)\mathcal{R}^{\epsilon}_{3}D^{(1)}_1D^{(1)}_1D^{(2)}_1P_0-f_{2, 1}(t, T)\mathcal{R}^{\epsilon}_{4}D^{(1)}_2D^{(2)}_1P_0\nonumber\\ & \quad -f_{0, 3}(t, T)\mathcal{R}^{\epsilon}_{5}D^{(2)}_1D^{(2)}_2P_0-f_{1, 2}(t, T)\mathcal{R}^{\epsilon}_{6}D^{(1)}_1D^{(2)}_1D^{(2)}_1P_0\nonumber\\ & \quad -f_{3, 0}(t, T)\mathcal{R}^{\epsilon}_{1}D^{(1)}_2P_0-f_{0, 3}(t, T)\mathcal{R}^{\epsilon}_{5}D^{(2)}_2P_0\nonumber\\ & \quad -f_{1, 2}(t, T)\mathcal{R}^{\epsilon}_{2}D^{(1)}_1D^{(2)}_1P_0-f_{2, 1}(t, T)\mathcal{R}^{\epsilon}_{4}D^{(1)}_1D^{(2)}_1P_0 \end{align}
with the terminal condition 
$\tilde {P}_{10}(T_0, x_1, x_2, z,T)=0$, where
\begin{equation} \begin{aligned} \mathcal{R}^{\epsilon}_1(z) & ={-}\frac{\sqrt{\epsilon}}{2}\rho_{1y}\zeta_1^3(z)\langle \zeta\beta \chi' \rangle,\\ \mathcal{R}^{\epsilon}_2(z) & ={-}\frac{\sqrt{\epsilon}}{2}\rho_{1y}\zeta_1(z)\zeta_2^2(z)\langle \zeta\beta \chi' \rangle,\\ \mathcal{R}^{\epsilon}_3(z) & ={-}\sqrt{\epsilon}\rho_{12}\rho_{1y}\zeta_1^2(z)\zeta_2(z)\langle \zeta\beta \chi' \rangle,\\ \mathcal{R}^{\epsilon}_4(z) & ={-}\frac{\sqrt{\epsilon}}{2}\rho_{2y}\zeta_1^2(z)\zeta_2(z)\langle \zeta\beta \chi' \rangle,\\ \mathcal{R}^{\epsilon}_5(z) & ={-}\frac{\sqrt{\epsilon}}{2}\rho_{2y}\zeta_2^3(z)\langle \zeta\beta \chi' \rangle,\\ \mathcal{R}^{\epsilon}_6(z) & ={-}\sqrt{\epsilon}\rho_{12}\rho_{2y}\zeta_1(z)\zeta_2^2(z)\langle \zeta\beta \chi' \rangle, \end{aligned} \end{equation}
\begin{align} & f_{i, j}(t, T)=e^{-(i\kappa_1+j\kappa_2)(T-t)}, \end{align}
\begin{align} & D^{(1)}_n=x_1^n\frac{\partial^n}{\partial x_1^n},\quad D^{(2)}_n=x_2^n\frac{\partial^n}{\partial x_2^n}, \end{align}
and 
$\chi =\chi (y)$ is the solution of the Poisson equation
\begin{equation} \mathcal{L}_0\chi =\zeta^2(y)-\langle \zeta^2 \rangle. \end{equation}
Then from the commutative relation between 
$D^{(j)}_n$ and 
$\mathcal {L}_{\mathrm {bln}}$, the solution 
$\tilde {P}_{10}$ is given by (42).
A.3. Derivation of solution $P_{01}$
From the results for 
$P_0$ and 
$P_{10}$ and the PDE 
$\mathcal {L}_0P_{21}+\mathcal {L}_1P_{11}+\mathcal {L}_2P_{01}+\mathcal {A}_3P_0+\mathcal {M}_1P_{10}=0$, one can find that 
$\tilde {P}_{01}=\sqrt {\delta }P_{01}$ satisfies
\begin{align} \langle \mathcal{L}_2 \rangle\tilde{P}_{01}& ={-}\sqrt{\delta}\langle \mathcal{A}_3 \rangle P_0\nonumber\\ & ={-}f_{1, 0}\mathcal{R}_1^\delta(z)D^{(1)}_1\frac{\partial P_0}{\partial z} - f_{0, 1}\mathcal{R}_2^\delta(z)D^{(2)}_1\frac{\partial P_0}{\partial z}\nonumber\\ & \quad -(f_{1, 0}-f_{3,0})\mathcal{R}_3^\delta(z)D^{(1)}_2P_0 - (f_{0, 1}-f_{0, 3})\mathcal{R}_4^\delta(z)D^{(2)}_2P_0\nonumber\\ & \quad -(f_{1, 0}-f_{1, 2})\mathcal{R}_5^\delta(z)D^{(1)}_1D^{(2)}_1P_0 - (f_{0, 1}-f_{2, 1})\mathcal{R}_6^\delta(z)D^{(1)}_1D^{(2)}_1P_0 \end{align}
with terminal condition 
$\tilde {P}_{01}(T_0, x_1, x_2, z,T)=0$, where 
$\mathcal {R}_i^\delta (z)$'s are given by
\begin{equation} \begin{aligned} & \mathcal{R}_1^\delta(z)=\sqrt{\delta}\rho_{1z}\langle\zeta\rangle\zeta_1(z)g(z),\\ & \mathcal{R}_2^\delta(z)=\sqrt{\delta}\rho_{2z}\langle\zeta\rangle\zeta_2(z)g(z),\\ & \mathcal{R}_3^\delta(z)=\sqrt{\delta}\frac{\rho_{1z}}{\kappa_1}\bar{\zeta}^2\langle\zeta\rangle\zeta_1^2(z)\zeta_1^\prime(z)g(z),\\ & \mathcal{R}_4^\delta(z)=\sqrt{\delta}\frac{\rho_{2z}}{\kappa_2}\bar{\zeta}^2\langle\zeta\rangle\zeta_2^2(z)\zeta_2^\prime(z)g(z),\\ & \mathcal{R}_5^\delta(z)=\sqrt{\delta}\frac{\rho_{1z}}{2\kappa_2}\bar{\zeta}^2\langle\zeta\rangle\zeta_1(z)\zeta_2(z)\zeta_2^\prime(z)g(z),\\ & \mathcal{R}_6^\delta(z)=\sqrt{\delta}\frac{\rho_{2z}}{2\kappa_1}\bar{\zeta}^2\langle\zeta\rangle\zeta_1(z)\zeta_1^\prime(z)\zeta_2(z)g(z), \end{aligned} \end{equation}
respectively. To solve the PDE problem (A.11), we utilize the relation
\begin{equation} \begin{aligned} \sigma_1\frac{\partial P}{\partial \sigma_1} & =\sigma_1^2x_1^2(T_0-t)\frac{\partial^2 P}{\partial x_1^2} + \rho\sigma_1\sigma_2x_1x_2(T_0-t)\frac{\partial^2 P}{\partial x_1 \partial x_2},\\ \sigma_2\frac{\partial P}{\partial \sigma_2} & =\sigma_2^2x_2^2(T_0-t)\frac{\partial^2 P}{\partial x_2^2} + \rho\sigma_1\sigma_2x_1x_2(T_0-t)\frac{\partial^2 P}{\partial x_1 \partial x_2} \end{aligned} \end{equation}
that hold for 
$P=P_{\mathrm {KBS}}$. This is the (two-dimensional version of) Vega-Gamma relation for spread option prices in the bivariate log-normal model. Then using the linearity of the differential operators involved, one can easily obtain
\begin{align} \tilde{P}_{01}& =\mu_1\mathcal{R}_1^\delta(z)D^{(1)}_1\frac{\partial P_0}{\partial z}+\mu_2\mathcal{R}_2^\delta(z)D^{(2)}_1\frac{\partial P_0}{\partial z}\nonumber\\ & \quad +(\nu_{1,0}-\nu_{3,0})\mathcal{R}_3^\delta(z)D^{(1)}_2P_0 +(\nu_{0,1}-\nu_{0,3})\mathcal{R}_4^\delta(z)D^{(2)}_2P_0\nonumber\\ & \quad +(\nu_{1,0}-\nu_{1,2})\mathcal{R}_5^\delta(z)D^{(1)}_1D^{(2)}_1P_0 +(\nu_{0,1}-\nu_{2,1})\mathcal{R}_6^\delta(z)D^{(1)}_1D^{(2)}_1P_0, \end{align}
where
\begin{equation} \begin{aligned} \mu_1 & =\frac{(T_0-t)f_{1,0}-\nu_{1,0}}{\kappa_1(T_0-t)},\\ \mu_2 & =\frac{(T_0-t)f_{0,1}-\nu_{0,1}}{\kappa_2(T_0-t)}. \end{aligned} \end{equation}
Open access
