Monte Carlo Simulation for Trading Under a Lévy-Driven Mean-Reverting Framework

Abstract

We present a Monte Carlo approach to pairs trading on mean-reverting spreads modelled by Lévy-driven Ornstein-Uhlenbeck processes. Specifically, we focus on using a variance gamma driving process, an infinite activity pure jump process to allow for more flexible models of the price spread than is available in the classical model. However, this generalization comes at the cost of not having analytic formulas, so we apply Monte Carlo methods to determine optimal trading levels and develop a variance reduction technique using control variates. Within this framework, we numerically examine how the optimal trading strategies are affected by the parameters of the model. In addition, we extend our method to bivariate spreads modelled using a weak variance alpha-gamma driving process, and explore the effect of correlation on these trades.

Keywords:

• Pairs trading
• Monte Carlo simulation
• Lévy process
• Ornstein-Uhlenbeck process
• mean reversion
• variance gamma process

1. Introduction

Many empirical studies have demonstrated examples of mean-reverting spreads in different markets. They can be generated by taking long and short positions in stocks and ETFs (Avellaneda and Lee 2010; Gatev, Goetzmann, and Rouwenhorst 2006; Huck and Afawubo 2015), futures contracts (Brennan and Schwartz; Dai, Zhong, and Kwok), physical commodity and commodity ETFs (Kanamura, Rachev, and Fabozzi 2010; Leung and Li 2016), as well as cryptocurrencies (Leung and Nguyen 2019).

Another popular approach, the stochastic spread method captures the path behaviour of the spread through a stochastic process with mean reversion. Given a mean-reverting spread process, we consider the problem of pairs trading, that is, extracting trading signals from the model based on the idea that a trade should be entered when the spread is far from the mean, and exited when it returns to the mean, so that the trader profits based on the tendency of the process to mean revert. This is a form of statistical arbitrage. The construction of spreads and extraction of trading signals are derived from the analysis of the underlying model. For instance, Elliott, Van Der Hoek, and Malcolm (2005) propose a mean-reverting Gaussian Markov chain model to describe spread dynamics. The model's estimates are compared with observations of the spread to determine appropriate trading decisions. They solve an optimal double-stopping problem to analyse the timing of entry and exit subject to transaction costs.

Most prominently, Ornstein-Uhlenbeck (OU) processes have been used to model spreads for pairs trading. For classical OU processes, which are driven by Brownian motion, closed-form solutions to the pairs trading problem were obtained in Leung and Li (2015) and Lipton and López de Prado (2020).

Lévy processes have been widely applied in finance to model price processes with jumps (Cont and Tankov 2004; Schoutens 2003), which leads us to model price spreads with jumps. These Lévy models often provide a better fit than models based on Brownian motion (Madan, Carr, and Chang 1998; Michaelsen and Szimayer 2018; Schoutens 2003), where the first two references focus on variance gamma-based models. Jumps can also arise from hard-to-borrow stocks (Avellaneda and Lipkin 2009) or illiquid stocks (Gardini, Sabino, and Sasso 2022). This motivates going beyond the classical OU processes to consider Lévy-driven Ornstein-Uhlenbeck processes (LDOUP) to model mean-reverting spreads with jumps in pairs trading.

LDOUPs have been applied to other areas of mathematical finance such as modelling stochastic volatility (Barndorff-Nielsen and Shephard 2001), energy prices and energy derivatives (Benth and Schmeck 2014; Cummins, Kiely, and Murphy 2018; Sabino 2020), credit risk (Cariboni and Schoutens 2009), among others. LDOUPs also have increased flexibility for capturing skewness, kurtosis, and dependence in the multivariate setting. Parametric estimation of LDOUPs using maximum likelihood has been considered in the univariate and multivariate cases by Valdivieso, Schoutens, and Tuerlinckx (2009) and Lu (2022), respectively, and these methods can be applied to fit LDOUPs to spreads. However, in the context of pairs trading, LDOUPs have been applied restricting the driving process to compound Poisson processes with double exponential jumps. This setup is used since Borovkov and Novikov (2008) derived an analytic formula for the expected exit time which relies on the memorylessness property of the exponential jumps to deal with overshoots. Without this assumption, overshoots present a difficult technical hurdle to finding analytic formulas for the passage times of these processes in general, which requires Monte Carlo methods to overcome. Based on Borovkov and Novikov (2008) and Wu, Zang, and Zhao (2020) provide an analytic value function for pairs trading, while Endres and Stübinger (2019) use the analytic formula of Borovkov and Novikov (2008), while ignoring overshoots, to conduct an empirical study of optimal pairs trading.

In this paper, we discuss a framework for pairs trading using Lévy-driven Ornstein-Uhlenbeck processes to model the mean-reverting spread, which works with any LDOUP that can be simulated. Specifically, we focus on a variance gamma driving process, which is a time-changed Brownian motion and a pure-jump, infinite activity Lévy process with skewness and kurtosis parameters. By using this driving process, these properties extend to the LDOUP, allowing for more flexible models of the price spread than is available in the classical OU model, and in contrast to the finite activity processes considered previously. However, the generalization comes at the cost of not having analytic formulas. Hence, we develop a Monte Carlo method to evaluate the trading strategies and determine the optimal trading levels. In particular, we introduce an array of control variates to reduce the variance of the Monte Carlo estimator of the value function. We derive a formula for the variance reduction ratio when the control variate is used to estimate the mean and variance of the profit simultaneously, which extends the usual approach where only one parameter is estimated by control variates. Furthermore, in our numerical implementation, we examine how the parameters of the model affect the optimal trading levels and optimal expected profit in various scenarios. Specifically, we demonstrate numerically how skewness in the distribution of the spread leads to asymmetry in the optimal trading levels, how the optimal exit level can differ from the stationary mean when the process starts sufficiently far away from it and the discount rate is large, and it is otherwise difficult to find other situations where this occurs. And lastly, we demonstrate how the jumps can increase the optimal entry level and optimal expected profit.

We note that all the studies mentioned above focus on the trading performance of a single spread. The recent paper by Lee, Leung, and Ning (2023) introduces a framework for trading multiple mean-reverting spreads simultaneously, where capital is dynamically allocated among different spreads based on their statistical characteristics. In this paper, we also consider the trading problem involving two mean-reverting spreads modelled by an LDOUP driven by a weak variance alpha-gamma process, which is a multivariate generalization of the variance gamma process introduced in Buchmann, Lu, and Madan (2019). For bivariate pairs trading, we estimate the expected profit via Monte Carlo simulation and examine the effect of the correlation between spreads on the optimal trading levels and optimal expected profit. We demonstrate an example where the optimal strategy is to effectively pick the best spread mostly ignoring correlation, and an example where having low absolute correlation is more profitable when monitoring multiple spreads with marginal components that follow the same law.

The rest of the paper is structured as follows. In Section 2, we review Lévy-driven Ornstein-Uhlenbeck processes and required properties used in this study. The trading problem is formulated in Section 3. We discuss our Monte Carlo estimation framework in Section 4 and present the numerical results in Section 5. In Section 6, we examine the trading problem involving two mean-reverting spreads. Concluding remarks are provided in Section 7.

2. Lévy-Driven Ornstein-Uhlenbeck Process

Let $\mathbf{Z}\sim L^n(\mathit{\mu },\mathrm{\Sigma },\mathcal{Z})$ denote an n-dimensional Lévy process with characteristic triplet $(\mathit{\mu },\mathrm{\Sigma },\mathcal{Z})$, where $\mathit{\mu }\in {\mathbb{R}}^n$, $\mathrm{\Sigma }\in {\mathbb{R}}^{n\times n}$ is a covariance matrix and $\mathcal{Z}$ is an n-dimensional Lévy measure. See Bertoin (1996), Sato (1999), and Schoutens (2003) for references on Lévy processes.

A Lévy-driven Ornstein-Uhlenbeck process (LDOUP) $\mathbf{X}$ is defined by the stochastic differential equation (1) $\mathrm{d}\mathbf{X}(t)=-\lambda \mathbf{X}(t)\mathrm{d}t+\mathrm{d}\mathbf{Z}(\mathrm{\lambda t}),\mathbf{X}(0)={\mathbf{X}}_0,t\ge 0,$(1) where $\mathbf{Z}\sim L^n(\mathit{\mu },\mathrm{\Sigma },\mathcal{Z})$ is the background driving Lévy process (BDLP), $\lambda >0$, and ${\mathbf{X}}_0$ is a random vector independent of $\mathbf{Z}$. Note that there is no loss in generality in using the drift term $-\lambda \mathbf{X}(t)\mathrm{d}t$ instead of $\lambda (\overline{\mathit{\mu }}-\mathbf{X}(t))\mathrm{d}t$ for some $\overline{\mathit{\mu }}\in {\mathbb{R}}^n$, since the stationary mean of $\mathbf{X}$ is controlled by $\mathbb{E}[\mathbf{Z}(1)]$. Also, there is no loss in generality in using the BDLP $\mathbf{Z}\circ (\mathrm{\lambda I})$, where I is the identity function, instead of $\mathbf{Z}$, since any Lévy process $\tilde{\mathbf{Z}}\sim L^n(\tilde{\mathit{\mu }},\tilde{\mathrm{\Sigma }},\tilde{\mathcal{Z}})$ can be written in the form $\tilde{\mathbf{Z}}\stackrel{D}{=}\mathbf{Z}\circ (\mathrm{\lambda I})$, where $\mathbf{Z}\sim L^n(\tilde{\mathit{\mu }}/\lambda ,\tilde{\mathrm{\Sigma }}/\lambda ,\tilde{\mathcal{Z}}/\lambda )$. However, the former is more convenient since it leads to a stationary distribution that does not depend on λ. The exposition on LDOUP here follows Lu (2022), while some additional references for LDOUPs include Masuda (2004), Sato and Yamazato (1983), and Sato and Yamazato (1985).

The solution of (1(1) $\mathrm{d}\mathbf{X}(t)=-\lambda \mathbf{X}(t)\mathrm{d}t+\mathrm{d}\mathbf{Z}(\mathrm{\lambda t}),\mathbf{X}(0)={\mathbf{X}}_0,t\ge 0,$(1) ) is (2) $\mathbf{X}(t)=e^{-\mathrm{\lambda t}}\mathbf{X}(0)+e^{-\mathrm{\lambda t}}{\int }_0^t{e}^{\mathrm{\lambda s}}\mathrm{d}\mathbf{Z}(\mathrm{\lambda s}),t\ge 0.$(2) For $t_0=0,t_1=\mathrm{\Delta },\dots ,t_q=q\mathrm{\Delta }$, at the observation time $t_i$, we have $\mathbf{X}(t_i)=e^{-\lambda \mathrm{\Delta }}(\mathbf{X}(t_{i-1})+{\int }_{t_{i-1}}^{t_i}{e}^{\mathrm{\lambda s}}\mathrm{d}\mathbf{Z}(\mathrm{\lambda s})),i=1,\dots ,q.$ The stochastic integral term is iid, and we let ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })$ be the random vector with this distribution, so it is interpreted as the innovation term, up to a factor $e^{-\lambda \mathrm{\Delta }}$, of an AR(1) process. Specifically, we define ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })={\int }_0^{\mathrm{\Delta }}{e}^{\mathrm{\lambda s}}\mathrm{d}\mathbf{Z}(\mathrm{\lambda s})\stackrel{D}{=}{\int }_{t_{i-1}}^{t_i}{e}^{\mathrm{\lambda s}}\mathrm{d}\mathbf{Z}(\mathrm{\lambda s}),\mathrm{\Delta }>0.$ The random vector ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })$ is infinitely divisible, so it is determined by its characteristic exponent, which is given by (3) ${\mathrm{\Psi }}_{{\mathbf{Z}}^{\ast }(\mathrm{\Delta })}(\mathit{\theta })={\int }_0^{\lambda \mathrm{\Delta }}{\mathrm{\Psi }}_{\mathbf{Z}}(e^t\mathit{\theta })\mathrm{d}t,\mathit{\theta }\in {\mathbb{R}}^n,$(3) where ${\mathrm{\Psi }}_{\mathbf{Z}}$ is the characteristic exponent of the BLDP $\mathbf{Z}$. Combining these results, the LDOUP can be written as (4) $\mathbf{X}(t_i)=e^{-\lambda \mathrm{\Delta }}(\mathbf{X}(t_{i-1})+{\mathbf{Z}}^{\ast }(\mathrm{\Delta }{)}^{(i)}),i=1,\dots ,q,$(4) where ${\mathbf{Z}}^{\ast }(\mathrm{\Delta }{)}^{(i)}$ is iid with the distribution of ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })$.

Next, we consider specific examples of LDOUP processes where the BDLP is a tractable pure-jump infinite activity Lévy process to model price spreads for pairs trading. We consider both the univariate and multivariate cases, where the BDLP is a variance gamma and a weak variance alpha-gamma process, respectively.

2.1. Univariate LDOUP Driven by a Variance Gamma Process

The variance gamma process was introduced in Madan and Seneta (1990) and additional details can be found in Madan, Carr, and Chang (1998). Let ${\mathrm{\Gamma }}_S(a,b)$ denote a gamma subordinator with shape a>0 and rate b>0, and $B{M}^n(\mathit{\mu },\mathrm{\Sigma })$ denote a Brownian motion with drift $\mathit{\mu }=({\mu }_1,\dots ,{\mu }_n)\in {\mathbb{R}}^n$ and covariance matrix $\mathrm{\Sigma }=({\mathrm{\Sigma }}_{\mathrm{kl}})\in {\mathbb{R}}^{n\times n}$. A Lévy process $\mathbf{V}\sim V{G}^n(b,\mathit{\mu },\mathrm{\Sigma })$ is a variance gamma (VG) process if $\mathbf{V}\stackrel{D}{=}\mathit{\eta I}+\mathbf{B}\circ (G,\dots ,G),$ where $\mathbf{B}\sim B{M}^n(\mathit{\mu },\mathrm{\Sigma })$ and $G\sim {\mathrm{\Gamma }}_S(b,b)$ are independent, $\mathit{\eta }\in {\mathbb{R}}^n$, and I is the identity function. The parameter $\mathit{\eta }$ adds a drift.

We call the univariate LDOUP $X\sim \mathrm{OU}\text{-}\mathrm{VG}(\lambda ,b,\mu ,{\sigma }^2,\eta )$ an $\mathrm{OU}\text{-}\mathrm{VG}$ process if its BDLP is $Z\sim V{G}^1(b,\mu ,{\sigma }^2,\eta )$.

2.2. Multivariate LDOUP Driven by a Weak Variance Alpha-Gamma Process

The weak variance alpha-gamma (WVAG) process was introduced in Buchmann, Lu, and Madan (2019). Let $n\ge 2$, a>0, $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_n)\in (0,1/a{)}^n$, ${\beta }_k:=(1-a{\alpha }_k)/{\alpha }_k$, $k=1,\dots ,n$ and $\mathit{\eta }\in {\mathbb{R}}^n$. Let $\mathit{\alpha }\diamond \mathit{\mu }:=({\alpha }_1{\mu }_1,\dots ,{\alpha }_n{\mu }_n)\in {\mathbb{R}}^n$ and $\mathit{\alpha }\diamond \mathrm{\Sigma }:=({\mathrm{\Sigma }}_{\mathrm{kl}}({\alpha }_k\wedge {\alpha }_l))\in {\mathbb{R}}^{n\times n}$. We define $\mathbf{Z}\sim \mathrm{WVA}{G}^n(a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$ as a weak variance alpha-gamma (WVAG) process if (5) $\mathbf{Z}\stackrel{D}{=}\mathit{\eta I}+{\mathbf{V}}_0+(V_1,\dots ,V_n),$(5) where ${\mathbf{V}}_0\sim V{G}^n(a,a\mathit{\alpha }\diamond \mathit{\mu },a\mathit{\alpha }\diamond \mathrm{\Sigma })$ and $V_k\sim V{G}^1({\beta }_k,{\alpha }_k{\beta }_k{\mu }_k,{\alpha }_k{\beta }_k{\mathrm{\Sigma }}_{\mathrm{kk}})$, $k=1,\dots ,n$, are independent. This formulation gives a direct method of simulating the WVAG process in terms of VG processes, and seeing that the process has common and idiosyncratic jumps.

The original definition of the WVAG process (with $\mathit{\eta }=\mathbf{0}$), (6) $\mathbf{Z}\stackrel{D}{=}\mathbf{B}\odot \mathbf{T},$(6) where $\mathbf{B}\sim B{M}^n(\mathit{\mu },\mathrm{\Sigma })$, $\mathbf{T}$ an n-dimensional alpha-gamma subordinator, and ⊙ is the weak subordination operation, gives a multivariate time-change interpretation to the process, generalizing the univariate VG process. Indeed, the marginal components of the WVAG process are general univariate VG processes with $Z_k\sim V{G}^1(1/{\alpha }_k,{\mu }_k,{\mathrm{\Sigma }}_{\mathrm{kk}},{\eta }_k)$.

We call a multivariate LDOUP $\mathbf{X}\sim \mathrm{OU}\text{-}\mathrm{WVA}{G}^n(\lambda ,a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$ an OU-WVAG process if its BDLP is $\mathbf{Z}\sim \mathrm{WVA}{G}^n(a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$.

Based on the moment formulas (see Michaelsen and Szimayer 2018), the parameters of the WVAG process can be interpreted as follows. The marginal parameters ${\alpha }_k,{\mu }_k,{\mathrm{\Sigma }}_{\mathrm{kk}},{\eta }_k$ are kurtosis, skewness, variance, and location parameters, respectively, that affect the marginal component distribution, while $a,{\mathrm{\Sigma }}_{\mathrm{kl}},k\ne l$, are joint parameters. In particular, as ${\alpha }_k\to 0$ (equivalently, the VG parameter $b_k=1/{\alpha }_k\to \mathrm{\infty }$ in the kth component), the kth marginal component of the WVAG process converges in law to Brownian motion so the sample paths become closer to continuous.

2.3. Moment Formulas

For a random vector $\mathbf{U}=(U_1,\dots ,U_n)$, denote the mean by $m_1(\mathbf{U}):=\mathbb{E}[\mathbf{U}]$, and the kth central moment of $\mathbf{U}$ by $m_k(\mathbf{U}):=(\mathbb{E}[(U_1-\mathbb{E}[U_1]{)}^k],\dots ,\mathbb{E}[(U_n-\mathbb{E}[U_n]{)}^k])$, $k\ge 2$.

For a general BDLP $\mathbf{Z}$, Lu (2022, Lemma 1) gives formulas for the moments of ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })$ in terms of the moments of $\mathbf{Z}$. Combining this with (2(2) $\mathbf{X}(t)=e^{-\mathrm{\lambda t}}\mathbf{X}(0)+e^{-\mathrm{\lambda t}}{\int }_0^t{e}^{\mathrm{\lambda s}}\mathrm{d}\mathbf{Z}(\mathrm{\lambda s}),t\ge 0.$(2) ) and conditional on ${\mathbf{X}}_0$ (or assuming it is constant), the moments of the LDOUP $\mathbf{X}$ are given by $\begin{array}{rl}{m}_1(\mathbf{X}(t))& =e^{-\mathrm{\lambda t}}{\mathbf{X}}_0+(1-e^{-\mathrm{\lambda t}})m_1(\mathbf{Z}(1)),\\ m_2(\mathbf{X}(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)m_2(\mathbf{Z}(1)),\\ m_3(\mathbf{X}(t))& =\left(\frac{1-e^{-3\mathrm{\lambda t}}}{3}\right)m_3(\mathbf{Z}(1)),\\ \mathrm{Cov}(X_k(t),X_l(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)\mathrm{Cov}(Z_k(1),Z_l(1)),k\ne l.\end{array}$ The stationary distribution $\mathbf{Y}$ of $\mathbf{X}$ is the distribution such that if ${\mathbf{X}}_0\stackrel{D}{=}\mathbf{Y}$, then $\mathbf{X}(t)\stackrel{D}{=}\mathbf{Y}$ for all $t\ge 0$. Thus, the stationary mean, which is the mean of the stationary distribution, is $\overline{\mathit{\mu }}=\mathbb{E}[\mathbf{Z}(1)]$. This is the level that we expect the process to return to in pairs trading.

Now specializing to the case where $\mathbf{X}\sim \mathrm{OU}\text{-}\mathrm{WVA}{G}^n(\lambda ,a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$, and using the moment formulas for the WVAG process in Michaelsen and Szimayer (2018, Remark 4 and Appdendix A.1), we have (7) $\begin{array}{rl}\mathbb{E}[X_1(t)]& =e^{-\mathrm{\lambda t}}{X}_0+(1-e^{-\mathrm{\lambda t}})({\eta }_1+{\mu }_1),\end{array}$(7) (8) $\begin{array}{rl}\mathrm{Var}(X_1(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)({\mathrm{\Sigma }}_{11}+{\alpha }_1{\mu }_1^2),\end{array}$(8) (9) $\begin{array}{rl}\mathrm{Skew}(X_1(t))& =\frac{(\frac{1-e^{-3\mathrm{\lambda t}}}{3})(3{\mathrm{\Sigma }}_{11}{\mu }_1{\alpha }_1+2{\mu }_1^3{\alpha }_1^2)}{{((\frac{1-e^{-2\mathrm{\lambda t}}}{2})({\mathrm{\Sigma }}_{11}+{\mu }_1^2{\alpha }_1))}^{3/2}}\end{array}$(9) (10) $\begin{array}{rl}\mathrm{Cov}(X_1(t),X_2(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)a(({\alpha }_1\wedge {\alpha }_2){\mathrm{\Sigma }}_{12}+{\alpha }_1{\alpha }_2{\mu }_1{\mu }_2).\end{array}$(10) Furthermore, $\overline{\mathit{\mu }}=\mathit{\mu }+\mathit{\eta }$. These formulas will be used for control variates in the Monte Carlo method. Note that (7(7) $\begin{array}{rl}\mathbb{E}[X_1(t)]& =e^{-\mathrm{\lambda t}}{X}_0+(1-e^{-\mathrm{\lambda t}})({\eta }_1+{\mu }_1),\end{array}$(7) ) to (9(9) $\begin{array}{rl}\mathrm{Skew}(X_1(t))& =\frac{(\frac{1-e^{-3\mathrm{\lambda t}}}{3})(3{\mathrm{\Sigma }}_{11}{\mu }_1{\alpha }_1+2{\mu }_1^3{\alpha }_1^2)}{{((\frac{1-e^{-2\mathrm{\lambda t}}}{2})({\mathrm{\Sigma }}_{11}+{\mu }_1^2{\alpha }_1))}^{3/2}}\end{array}$(9) ) apply to OU-VG processes with $b_1=1/{\alpha }_1$.

3. Pairs Trading Problem

Suppose the spread of a pair of assets or the price of a portfolio of assets is mean-reverting and follows a univariate LDOUP X.

Now we proceed to determine the levels at which trades should be entered or exited. Consider the following pairs trading strategy on the spread X: Enter a short position in X when it is above $\overline{\mu }+d$, and close the position by buying when it is below $\overline{\mu }+c$, where $0\le c<d$. Enter a long position when it is below $\overline{\mu }-d$, and close the position by selling when it is above $\overline{\mu }-c$. Fix a terminal time T>0, where the trade is closed by time T if it has not been closed already. Thus, d and c are the entry and exit levels relative to the mean reversion level $\overline{\mu }$.

This trading strategy is depicted for a sample path of X in Figure 1. Note that the price at which the trade enters and exits is not the same as the entry and exit levels due to the presence of jumps. In this model when X is an OU-VG process, there is always some degree of overshoot.

Figure 1. For the spread process X (black line), a trade is entered when the price passes $\pm d$, whichever happens first (red lines), and exited when it passes $\pm c=0$ (blue line). The value of the spread at the entry (red point) and exit (blue point) times are shown. Here, $X\sim \mathrm{OU}\text{-}\mathrm{VG}(1,5,0,0.015,0)$.

More generally, we can consider possibly different entry levels $d_{+}$, $d_{-}$ from above and below $\overline{\mu }$, and their respective exit levels $c_{+},c_{-}$. These satisfy $0\le c_{\pm }<d_{\pm }$.

The relevant passage times are $\begin{array}{rl}{\tau }_{d_{+}}& =inf\{t\ge 0:X(t)>\overline{\mu }+d_{+}\}\wedge T,\\ {\tau }_{c_{+}}& =inf\{t\ge {\tau }_{d_{+}}:X(t)<\overline{\mu }+c_{+}\}\wedge T,\\ {\tau }_{d_{-}}& =inf\{t\ge 0:X(t)<\overline{\mu }-d_{-}\}\wedge T,\\ {\tau }_{c_{-}}& =inf\{t\ge {\tau }_{d_{-}}:X(t)>\overline{\mu }-c_{-}\}\wedge T.\end{array}$ Let $r\in \mathbb{R}$ be a discount rate, then P is the discounted profit over one trade cycle given by (11) $\begin{array}{rl}P:=P(c_{+},c_{-},d_{+},d_{-})& ={\mathbf{1}}_{\{{\tau }_{d_{+}}\le {\tau }_{d_{-}}\}}{e}^{-r{\tau }_{c_{+}}}(X({\tau }_{d_{+}})-X({\tau }_{c_{+}}))\\ & +{\mathbf{1}}_{\{{\tau }_{d_{+}}>{\tau }_{d_{-}}\}}{e}^{-r{\tau }_{c_{-}}}(X({\tau }_{c_{-}})-X({\tau }_{d_{-}})).\end{array}$(11) For convenience, we will refer to this simply as the profit. The value function is (12) $\begin{array}{rl}V:=V(c_{+},c_{-},d_{+},d_{-})& =\mathbb{E}[P]-\mathrm{\gamma Var}(P)\\ & =\mathbb{E}[P]-\gamma \mathbb{E}[P^2]+\gamma \mathbb{E}[P{]}^2,\end{array}$(12) where $\gamma \ge 0$ is a parameter which penalizes trading strategies with high variance. If $\gamma =0$, the value function is the expected profit.

We want to find the optimal levels $c_{\pm },d_{\pm }$ for the pairs trading strategy to maximize the value function, that is, $\underset{c_{+},c_{-},d_{+},d_{-}}{max}V(c_{+},c_{-},d_{+},d_{-}).$ We develop a Monte Carlo approach to simulate sample paths of the LDOUP X, to estimate $V(c_{+},c_{-},d_{+},d_{-})$ on $0\le c_{\pm }<d_{\pm }$, and in turn numerically solve the maximization problem for the optimal trading levels. Given the use of simulation, we also discuss control variates for variance reduction.

4. Monte Carlo Method

Let m be the number of Monte Carlo simulations. For a spread process following a univariate LDOUP X, suppose we have m simulations of $X(t_0),\dots ,X(t_q)$ with a small stepsize $\mathrm{\Delta }>0$. Then we can compute the random variable P, and hence obtain the Monte Carlo estimate of $\mathbb{E}[P]$ and $\mathbb{E}[P^2]$ to estimate the value function V in (12(12) $\begin{array}{rl}V:=V(c_{+},c_{-},d_{+},d_{-})& =\mathbb{E}[P]-\mathrm{\gamma Var}(P)\\ & =\mathbb{E}[P]-\gamma \mathbb{E}[P^2]+\gamma \mathbb{E}[P{]}^2,\end{array}$(12) ). The same sample paths are used for all evaluations of the same value function. While the Monte Carlo method presented here works for any LDOUP that can be fitted and simulated, we specifically focus on the OU-VG process.

4.1. Simulating the Spread

Let the spread process be $X\sim \mathrm{OU}\text{-}\mathrm{VG}(\lambda ,b,\mu ,{\sigma }^2,\eta )$. We assume that the parameters are known. For example, they may be estimated using maximum likelihood as in Valdivieso, Schoutens, and Tuerlinckx (2009) and Lu (2022).

It is possible to exactly simulate $Z^{\ast }(\mathrm{\Delta })$ for the OU-VG process, and hence $X(t_0),,\dots ,X(t_q)$ by (4(4) $\mathbf{X}(t_i)=e^{-\lambda \mathrm{\Delta }}(\mathbf{X}(t_{i-1})+{\mathbf{Z}}^{\ast }(\mathrm{\Delta }{)}^{(i)}),i=1,\dots ,q,$(4) ), using the method of Sabino (2020) and Qu, Dassios, and Zhao (2021). For $G\sim {\mathrm{\Gamma }}_S(a,b)$, define the random variable $G^{\ast }(\mathrm{\Delta })\stackrel{D}{=}{e}^{\lambda \mathrm{\Delta }}(\tilde{G}+\tilde{C})$, where $\tilde{G}\sim \mathrm{\Gamma }(\mathrm{a\lambda }\mathrm{\Delta },b{e}^{\lambda \mathrm{\Delta }})$ is gamma where the parameters are the shape and rate, respectively, and $\tilde{C}\sim \mathrm{CP}(a{\lambda }^2{\mathrm{\Delta }}^2/2,{\mathcal{P}}_J)$ is compound Poisson where the parameters are the rate and jump distribution, respectively. Here, ${\mathcal{P}}_J$ is the probability law of J such that $J|U\sim \mathrm{Exponential}(b{e}^{\lambda \mathrm{\Delta }\sqrt{U}})$ where the parameter is the rate, and $U\sim \mathrm{Uniform}(0,1)$. Then $Z^{\ast }(\mathrm{\Delta })\stackrel{D}{=}\eta (e^{\lambda \mathrm{\Delta }}-1)+G_{+}^{\ast }(\mathrm{\Delta })-G_{-}^{\ast }(\mathrm{\Delta }),$ where $G_{+}\sim {\mathrm{\Gamma }}_S(b,b_{+})$, $G_{-}\sim {\mathrm{\Gamma }}_S(b,b_{-})$, and $b_{\pm }:=\frac{2b}{\sqrt{{\mu }^2+2{\sigma }^{2}b}\pm \mu }.$

4.2. Control Variates

Next, we use control variates to reduce the variance of the Monte Carlo estimators of the value function. We refer to Glasserman (2003, Chapter 4) for the details of the method. Specifically, to estimate $\mathbb{E}[P]$, the method of control variates defines the random variable $P_C=P-\beta (C-\mathbb{E}[C])$, where C is the control variate for which $\mathbb{E}[C]$ has an analytic representation, and β is chosen to minimize $\mathrm{Var}(P_C)$, resulting in an unbiased estimator with reduced variance. Since the true value of β is unknown, the control variate estimator ${\hat{P}}_C$ is the sample mean of $P_C$ with β estimated by $\hat{\beta }$, the least squares estimator (LSE) of the slope in a linear regression where C is the regressor and P is the response variable. Note that $\hat{\beta }$ is also the sample estimate of the true β.

This method can be extended to multiple control variates, and the control variate estimator is equivalent to the prediction of a linear regression. See Appendix for details.

We apply the control variate method to estimate both $\mathbb{E}[P]$ and $\mathbb{E}[P^2]$ in the value function (12(12) $\begin{array}{rl}V:=V(c_{+},c_{-},d_{+},d_{-})& =\mathbb{E}[P]-\mathrm{\gamma Var}(P)\\ & =\mathbb{E}[P]-\gamma \mathbb{E}[P^2]+\gamma \mathbb{E}[P{]}^2,\end{array}$(12) ). Consider the control variates $X(t_1)-\overline{\mu },\dots ,X(t_p)-\overline{\mu }$, where $0,t_1,\dots ,t_p$ are equally spaced points on $[0,T]$. When estimating $\mathbb{E}[P^2]$, we include the additional control variates $(X(t_1)-\overline{\mu }{)}^2,\dots ,(X(t_p)-\overline{\mu }{)}^2$. Note that $\mathbb{E}[X(t)-\overline{\mu }]$ and $\mathbb{E}[(X(t)-\overline{\mu }{)}^2]$ can be analytically computed using (7(7) $\begin{array}{rl}\mathbb{E}[X_1(t)]& =e^{-\mathrm{\lambda t}}{X}_0+(1-e^{-\mathrm{\lambda t}})({\eta }_1+{\mu }_1),\end{array}$(7) ) and (8(8) $\begin{array}{rl}\mathrm{Var}(X_1(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)({\mathrm{\Sigma }}_{11}+{\alpha }_1{\mu }_1^2),\end{array}$(8) ).

Further, we may want to consider the event of entering and exiting a trade, and the event of entering and not exiting over the observations $X(t_1),\dots ,X(t_p)$ as additional control variates to approximate the indicator variable appearing in the profit function (11(11) $\begin{array}{rl}P:=P(c_{+},c_{-},d_{+},d_{-})& ={\mathbf{1}}_{\{{\tau }_{d_{+}}\le {\tau }_{d_{-}}\}}{e}^{-r{\tau }_{c_{+}}}(X({\tau }_{d_{+}})-X({\tau }_{c_{+}}))\\ & +{\mathbf{1}}_{\{{\tau }_{d_{+}}>{\tau }_{d_{-}}\}}{e}^{-r{\tau }_{c_{-}}}(X({\tau }_{c_{-}})-X({\tau }_{d_{-}})).\end{array}$(11) ). But the event that $X(t_1),\dots ,X(t_p)$ passes some level cannot be computed analytically for use as a control variate. Instead, the probability of a related event using the iid sequence of innovations $Z^{\ast }(\mathrm{\Delta }{)}^{(i)}$, $i=1,\dots ,p$ can be. Specifically, this motivates the inclusion of two additional control variates ${\mathbf{1}}_A$ and ${\mathbf{1}}_B$. Define A as the event that $Z^{\ast }(\mathrm{\Delta }{)}^{(1)},\dots ,Z^{\ast }(\mathrm{\Delta }{)}^{(p)}$ “enters a trade” (even though there is no actual trading on this sequence) by passing the level $e^{\lambda \mathrm{\Delta }}(\overline{\mu }+d)-\overline{\mu }$ or $e^{\lambda \mathrm{\Delta }}(\overline{\mu }-d)-\overline{\mu }$ and then “exits the trade” by passing the level $e^{\lambda \mathrm{\Delta }}(\overline{\mu }+c)-\overline{\mu }$ or $e^{\lambda \mathrm{\Delta }}(\overline{\mu }-c)-\overline{\mu }$, respectively.

The modified levels for $Z^{\ast }(\mathrm{\Delta }{)}^{(i)}$ are the levels such that $X(t_i)=e^{-\lambda \mathrm{\Delta }}(X(t_{i-1})+Z^{\ast }(\mathrm{\Delta }{)}^{(i)})$ passes the entry levels $\mu \pm d$ when $X(t_{i-1})=\overline{\mu }$, and similarly for the exit levels. Further, B is the same as A except without exit. Thus, A and B are the events that $Z^{\ast }(\mathrm{\Delta }{)}^{(i)}$, $i=1,\dots ,p$, enters and exits these modified levels, and enters and does not exit, respectively. Due to collinearity, the complementary event that the sequence does not enter is not included, and if p = 1, the event A, which has probability 0, is also not included. We have expressed this for symmetric levels $c=c_{+}=c_{-}$ and $d=d^{+}=d_{-}$, but this can be modified in an obvious way for the asymmetric case.

Since $Z^{\ast }(\mathrm{\Delta }{)}^{(i)}$, $i=1,\dots ,p$, are iid, we can compute the probabilities of A and B since we can compute $\mathbb{P}(Z^{\ast }(\mathrm{\Delta })\le x)$ by applying the Fourier inversion to the characteristic function of $Z^{\ast }(\mathrm{\Delta })$ determined by (3(3) ${\mathrm{\Psi }}_{{\mathbf{Z}}^{\ast }(\mathrm{\Delta })}(\mathit{\theta })={\int }_0^{\lambda \mathrm{\Delta }}{\mathrm{\Psi }}_{\mathbf{Z}}(e^t\mathit{\theta })\mathrm{d}t,\mathit{\theta }\in {\mathbb{R}}^n,$(3) ). The Fourier inversion methodology is described in Michaelsen and Szimayer (2018, Section 4.1).

4.3. Variance Reduction

Consider the two separate linear regression models, (13) ${\mathbf{Y}}_1=X_1{\mathit{\beta }}_1+{\mathit{ϵ}}_1,{\mathbf{Y}}_2=X_2{\mathit{\beta }}_2+{\mathit{ϵ}}_2,$(13) where ${\mathbf{Y}}_1$ and ${\mathbf{Y}}_2$ are the vector of m simulated values of P and $P^2$, respectively. In this section only, we let $X_1\in {\mathbb{R}}^{(p_1+1)\times m}$ and $X_2\in {\mathbb{R}}^{(p_2+1)\times m}$ be the design matrices, including an intercept. As outlined in Section 4.2, $X_1$ includes the control variates $X(t_i)-\overline{\mu }$, $i=1,\dots ,p$, ${\mathbf{1}}_A$, ${\mathbf{1}}_B$, and $X_2$ includes the control variates in $X_1$ in addition to $(X(t_i)-\overline{\mu }{)}^2$, $i=1,\dots ,p$. Since each simulated sample path of the LDOUP is independent, we have $\mathrm{Cov}({\mathit{ϵ}}_i)={\mathrm{\Sigma }}_{\mathit{ϵ},\mathrm{ii}}I$, i = 1, 2 and $\mathrm{Cov}({\mathit{ϵ}}_1,{\mathit{ϵ}}_2)={\mathrm{\Sigma }}_{\mathit{ϵ},12}I$ for some $2\times 2$ covariance matrix ${\mathrm{\Sigma }}_{\mathit{ϵ}}=({\mathrm{\Sigma }}_{\mathit{ϵ},\mathrm{ij}}{)}_{i,j=1}^2$, where I is the identity matrix.

For each of the linear regression models, i = 1, 2, let ${\mathit{\mu }}_{C,i}$ be the mean vector of the control variates computed as outlined in Section 4.2, then the predictions ${\mathbf{x}}_i=(1,{\mathit{\mu }}_{C,i})$ are ${\hat{Y}}_i={\mathbf{x}}_i^{\mathrm{\prime }}{\hat{\mathit{\beta }}}_i$. As explained in Appendix, the control variate estimators of $\mathbb{E}[P]$ and $\mathbb{E}[P^2]$ are the predictions ${\hat{Y}}_1$ and ${\hat{Y}}_2$, respectively. Thus, for each argument, the control variate estimator of the value function in (12(12) $\begin{array}{rl}V:=V(c_{+},c_{-},d_{+},d_{-})& =\mathbb{E}[P]-\mathrm{\gamma Var}(P)\\ & =\mathbb{E}[P]-\gamma \mathbb{E}[P^2]+\gamma \mathbb{E}[P{]}^2,\end{array}$(12) ) is ${\hat{V}}_C:={\hat{Y}}_1-\gamma {\hat{Y}}_2+\gamma {\hat{Y}}_1^2={\mathbf{b}}^{\mathrm{\prime }}\hat{\mathbf{Y}}+{\hat{\mathbf{Y}}}^{\mathrm{\prime }}A\hat{\mathbf{Y}},$ where $A=\left(\begin{array}{cc}\gamma & 0\\ 0& 0\end{array}\right),\mathbf{b}=\left(\begin{array}{c}1\\ -\gamma \end{array}\right),\hat{\mathbf{Y}}=\left(\begin{array}{c}{\hat{Y}}_1\\ {\hat{Y}}_2\end{array}\right).$ Let ${\overline{Y}}_1,{\overline{Y}}_2$ be the sample mean of ${\mathbf{Y}}_1,{\mathbf{Y}}_2$ respectively, then ${\hat{V}}_M:={\overline{Y}}_1-\gamma {\overline{Y}}_2+\gamma {\overline{Y}}_1^2$ is the corresponding Monte Carlo estimator of the value function.

Given that our approach involves parameters that are estimated by a function of two control variate estimators, rather than the usual situation of only one, we now derive the variance reduction factor in this context. For a large number of Monte Carlo simulations m, the asymptotic normality of the LSE gives $\hat{\mathbf{Y}}\sim N({\mathit{\mu }}_C,{\mathrm{\Sigma }}_C)$ approximately for some ${\mathit{\mu }}_C\in {\mathbb{R}}^2$ and ${\mathrm{\Sigma }}_C\in {\mathbb{R}}^{2\times 2}$. Using results on multivariate linear regression, and the law of total variance, noting that $X_1$ and $X_2$ are random, we have $\begin{array}{rl}{\mathit{\mu }}_C& =\left(\begin{array}{c}{\mathbf{x}}_1^{\mathrm{\prime }}{\mathit{\beta }}_1\\ {\mathbf{x}}_2^{\mathrm{\prime }}{\mathit{\beta }}_2\end{array}\right),\\ {\mathrm{\Sigma }}_C& =\left(\begin{array}{cc}{\mathrm{\Sigma }}_{\mathit{ϵ},11}{\mathbf{x}}_1^{\mathrm{\prime }}\mathbb{E}[(X_1^{\mathrm{\prime }}{X}_1{)}^{-1}]{\mathbf{x}}_1& {\mathrm{\Sigma }}_{\mathit{ϵ},12}{\mathbf{x}}_1^{\mathrm{\prime }}\mathbb{E}[(X_1^{\mathrm{\prime }}{X}_1{)}^{-1}{X}_1^{\mathrm{\prime }}{X}_2(X_2^{\mathrm{\prime }}{X}_2{)}^{-1}]{\mathbf{x}}_2\\ \ast & {\mathrm{\Sigma }}_{\mathit{ϵ},22}{\mathbf{x}}_2^{\mathrm{\prime }}\mathbb{E}[(X_2^{\mathrm{\prime }}{X}_2{)}^{-1}]{\mathbf{x}}_2\end{array}\right),\end{array}$ where $\ast$ denotes the rest of the matrix is completed by symmetry. By the plug-in principle, these parameters are estimated using $\begin{array}{rl}{\hat{\mathit{\mu }}}_C& =\left(\begin{array}{c}{\mathbf{x}}_1^{\mathrm{\prime }}{\hat{\mathit{\beta }}}_1\\ {\mathbf{x}}_2^{\mathrm{\prime }}{\hat{\mathit{\beta }}}_2\end{array}\right),\\ {\hat{\mathrm{\Sigma }}}_C& =\left(\begin{array}{cc}{\hat{\mathrm{\Sigma }}}_{\mathit{ϵ},11}{\mathbf{x}}_1^{\mathrm{\prime }}(X_1^{\mathrm{\prime }}{X}_1{)}^{-1}{\mathbf{x}}_1& {\hat{\mathrm{\Sigma }}}_{\mathit{ϵ},12}{\mathbf{x}}_1^{\mathrm{\prime }}(X_1^{\mathrm{\prime }}{X}_1{)}^{-1}{X}_1^{\mathrm{\prime }}{X}_2(X_2^{\mathrm{\prime }}{X}_2{)}^{-1}{\mathbf{x}}_2\\ \ast & {\hat{\mathrm{\Sigma }}}_{ϵ,22}{\mathbf{x}}_2^{\mathrm{\prime }}(X_2^{\mathrm{\prime }}{X}_2{)}^{-1}{\mathbf{x}}_2\end{array}\right),\end{array}$ where ${\mathrm{\Sigma }}_{\mathit{ϵ},\mathrm{ij}}$ is estimated using ${\hat{\mathrm{\Sigma }}}_{\mathit{ϵ},\mathrm{ij}}=\frac{{\widehat{\mathit{ϵ}}}_i^{\mathrm{\prime }}{\widehat{\mathit{ϵ}}}_j}{m-r-1},$ and ${\widehat{\mathit{ϵ}}}_1,{\widehat{\mathit{ϵ}}}_2$ are the residuals from the linear regressions in (13(13) ${\mathbf{Y}}_1=X_1{\mathit{\beta }}_1+{\mathit{ϵ}}_1,{\mathbf{Y}}_2=X_2{\mathit{\beta }}_2+{\mathit{ϵ}}_2,$(13) ), $r=p_1$ for i = j = 1, and $r=p_2$ otherwise.

Now since ${\hat{V}}_C$ is a quadratic form in $\hat{\mathbf{Y}}$, and $\hat{\mathbf{Y}}\sim N({\mathit{\mu }}_C,{\mathrm{\Sigma }}_C)$ holds approximately for large m, by Rencher and Schaalje (2008, Theorem 5.2c–d), the variance of the control variate estimator is (14) $\begin{array}{rl}\mathrm{Var}({\hat{V}}_C)& \approx {\mathbf{b}}^{\mathrm{\prime }}{\mathrm{\Sigma }}_C\mathbf{b}+2\mathrm{tr}((A{\mathrm{\Sigma }}_C{)}^2)+4{\mathit{\mu }}^{\mathrm{\prime }}A{\mathrm{\Sigma }}_{C}A{\mathit{\mu }}_C+4{\mathbf{b}}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{C}A{\mathit{\mu }}_C\\ & =(1+4{\mu }_{C,1}\gamma +4({\mu }_{C,1}{)}^2{\gamma }^2){\mathrm{\Sigma }}_{C,11}+2{\gamma }^2({\mathrm{\Sigma }}_{C,11}{)}^2\\ & -2\gamma (1+2{\mu }_{C,1}\gamma ){\mathrm{\Sigma }}_{C,12}+{\gamma }^2{\mathrm{\Sigma }}_{C,22}.\end{array}$(14) Then estimated variance $\widehat{\mathrm{Var}}({\hat{V}}_C)$ is obtained by replacing ${\mathit{\mu }}_C$ and ${\mathrm{\Sigma }}_C$ with ${\hat{\mathit{\mu }}}_C$ and ${\hat{\mathrm{\Sigma }}}_C$, respectively. By the asymptotic normality of ${\hat{V}}_M$, (14(14) $\begin{array}{rl}\mathrm{Var}({\hat{V}}_C)& \approx {\mathbf{b}}^{\mathrm{\prime }}{\mathrm{\Sigma }}_C\mathbf{b}+2\mathrm{tr}((A{\mathrm{\Sigma }}_C{)}^2)+4{\mathit{\mu }}^{\mathrm{\prime }}A{\mathrm{\Sigma }}_{C}A{\mathit{\mu }}_C+4{\mathbf{b}}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{C}A{\mathit{\mu }}_C\\ & =(1+4{\mu }_{C,1}\gamma +4({\mu }_{C,1}{)}^2{\gamma }^2){\mathrm{\Sigma }}_{C,11}+2{\gamma }^2({\mathrm{\Sigma }}_{C,11}{)}^2\\ & -2\gamma (1+2{\mu }_{C,1}\gamma ){\mathrm{\Sigma }}_{C,12}+{\gamma }^2{\mathrm{\Sigma }}_{C,22}.\end{array}$(14) ) is also used to compute $\widehat{\mathrm{Var}}({\hat{V}}_M)$ by replacing ${\hat{\mathit{\mu }}}_C$ and ${\hat{\mathrm{\Sigma }}}_C$ with ${\hat{\mathit{\mu }}}_M=({\overline{Y}}_1,{\overline{Y}}_2)$ and ${\hat{\mathrm{\Sigma }}}_M={\hat{\mathrm{\Sigma }}}_Y/m$, respectively, where ${\hat{\mathrm{\Sigma }}}_Y\in {\mathbb{R}}^{2\times 2}$ with the sample covariance $\widehat{\mathrm{Cov}}({\mathbf{Y}}_i,{\mathbf{Y}}_j)$ as its $(i,j)$-entry.

Thus, the estimated variance reduction factor of using control variates relative to Monte Carlo without control variates is $R=\frac{\widehat{\mathrm{Var}}({\hat{V}}_C)}{\widehat{\mathrm{Var}}({\hat{V}}_M)}.$ In the context of a single parameter estimated using control variates, while the true variance reduction factor using the true value of ${\mathit{\beta }}_i$ is no greater than 1, when using LSE ${\hat{\mathit{\beta }}}_i$ instead, R>1 is possible. This topic is discussed in Glasserman (2003, Section 4.1.3).

5. Numerical Implementation and Analysis

In this section, we implement the Monte Carlo method outlined above for various scenarios. The spread $X\sim \mathrm{OU}\text{-}\mathrm{VG}(\lambda ,b,\mu ,{\sigma }^2,\eta )$ is simulated with stepsize $\mathrm{\Delta }=0.01$ and terminal time T = 50. In all cases, $\mu =-\eta$, and any changes to μ adjusts η accordingly to ensure that the stationary mean is $\overline{\mu }=0$. Unless otherwise stated, the initial value is $X_0=0$, $\gamma =0$, the discount rate is r = 0.01, and it is assumed for simplicity that the entry level $d=d_{+}=d_{-}$ is symmetric and the exit level is fixed at $c=c_{+}=c_{-}=0$. Thus, we optimize for the entry level d. All results use m = 10, 000 Monte Carlo simulations.

5.1. Effect of Control Variates

Consider two examples using control variates.

Example 5.1

Suppose the spread X has parameters $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,1,-0.5,0.015,0.5)$, and $\gamma =0.1$. The number of time points p used in the control variate method can be chosen to obtain the optimal variance reduction factor, however, this also depends on the entry level d. We optimize this iteratively by using no control variates to optimize d, then use this value of d to optimize p with the resulting optimal value denoted $p^{\ast }$. Doing this, Figure 2 shows that for d = 1.086, the optimal number of time points for the control variates is $p^{\ast }=130$ when searching over $p=10,20,\dots ,200$. This gives a moderate optimal variance reduction factor of $R^{\ast }=0.735$.

Figure 2. Estimated variance reduction factor R as a function of the number of time points for the control variates p. The parameters are $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,1,-0.5,0.015,0.5)$, $X_0=0$, $\gamma =0.1$, r = 0.01.

Using $p^{\ast }=130$, Figure 3 shows the estimate of the value function $V(d)$, which has optimal entry level $d^{\ast }=1.086$. The control variate reduces $\mathrm{Var}({\hat{V}}_C)$, however, it does not change the optimal entry level $d^{\ast }$ in this example because the variability is already quite low, in fact, the coefficient of variation of ${\hat{V}}_C(d^{\ast })$, the estimate of $V(d^{\ast })$ conditional on $d^{\ast }$, is 0.0038 after control variates.

Figure 3. The estimate ${\hat{V}}_C(d)$ of the value function $V(d)$ and the optimal entry level $d^{\ast }=1.086$. The parameters are $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,1,-0.5,0.015,0.5)$, $X_0=0$, $\gamma =0.1$, r = 0.01.

Different choices in the parameters result in different variance reductions. In Table 1, results are shown varying $b,\mu$ (which also affects η), and for the parameters, decreasing b>0 increases $R^{\ast }$, while decreasing $\mu <0$ reduces and then increases $R^{\ast }$.

Table 1. For various values of $(b,\mu )$, the optimal number of time points for the control variates $p^{\ast }$ (searching over $p=10,20\dots ,200$) and the optimal variance reduction factor $R^{\ast }$ are shown.

b	μ	$p^{\ast }$	$R^{\ast }$
3	$-$0.05	70	0.951
2	$-$0.05	120	0.904
1	$-$0.05	100	0.871
1	$-$0.5	130	0.735
1	$-$1	110	0.826

Notes: The other parameters are $(\lambda ,{\sigma }^2)=(1,0.015)$, $X_0=0$, $\gamma =0.1$, r = 0.01.

Consider the case shown in the second row of Table 1 where $(b,\mu )=(2,-0.05)$ and $p^{\ast }=120$. Surprisingly, for large γ, such as $\gamma =1.5$, it is possible that the control variate method reduces the variance of the estimator of $\mathbb{E}[P]$ with an estimated variance reduction factor 0.910, and $\mathbb{E}[P^2]$ with estimated variance reduction factor 0.816, but increases the variance of the estimator of $V(d^{\ast })=\mathbb{E}[P]-\gamma \mathbb{E}[P^2]+\gamma \mathbb{E}[P{]}^2$ with variance reduction factor 1.139. In comparison, when $\gamma =0.1$, this does not occur, and those estimated variance reduction factors are 0.892, 0.793, 0.904, respectively.

Example 5.2

We use the same setup as Example 5.1, except X now has parameters $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(0.01,50,0.5,4,-0.5)$. Since $\lambda \approx 0$, this process has little mean reversion, and is close in law to a variance gamma process, so the value function is estimated with a relatively large amount of simulation error as shown in Figure 4. To address this, a locally estimated scatterplot smoothing (loess) smoother is applied to the estimated value function before maximizing it. Then the optimal entry level is $d^{\ast }=0.506$ without using control variates, and $d^{\ast }=0.456$ using the optimal number of control variates, which is $p^{\ast }=1$. The optimal variance reduction factor is $R^{\ast }=0.817$ and the coefficient of variation of ${\hat{V}}_C(d^{\ast })$ is 0.0281, which is consistent with the much larger variability of the estimated value function compared to Example 5.1, Figure 3. Unlike in Figure 3, using control variates here has an effect on the estimate of $d^{\ast }$, the size of this effect varies with the simulation. Note that $R^{\ast }$ is the variance reduction of estimating $V(d)$ for fixed d using control variates, and the method is not designed to reduce the standard error of $d^{\ast }$, although reducing the former should be expected to reduce the latter. For these optimal values, only 63.91% of simulations resulted in a trade that is entered and exited before the terminal time due to the slow mean reversion, compared to 90.88% in Figure 3.

Figure 4. The estimate of the value function $V(d)$ (solid lines) and the corresponding loess smooth (dashed lines) with (red lines) and without (black lines) control variates. The optimal entry level is $d^{\ast }=0.506$ without using control variates, and $d^{\ast }=0.456$ using the optimal number of control variates $p^{\ast }=1$. The parameters are $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(0.01,50,0.5,4,-0.5)$, $X_0=0$, $\gamma =0.1$, r = 0.01.

5.2. Asymmetry of Optimal Trading Levels

From now on, we set $\gamma =0$ and do not include control variates. In the previous examples, we fixed $d=d_{+}=d_{-}$ for simplicity. The parameter μ controls the asymmetry of the stationary distribution Y of the LDOUP X. If $\mu =0$, then the stationary distribution is symmetric and hence the optimal levels are equal, $d_{+}^{\ast }=d_{-}^{\ast }$. However, under asymmetry, we may want to optimize both $d_{+}$ and $d_{-}$ separately.

Here, we consider the optimal level in four cases: from a large negative skew to no skew. The other parameters of X are $(\lambda ,b,{\sigma }^2)=(1,5,0.015)$. The results are shown in Table 2, where $\mathrm{Skew}(Y)$ is the skewness of the stationary distribution, which by (9(9) $\begin{array}{rl}\mathrm{Skew}(X_1(t))& =\frac{(\frac{1-e^{-3\mathrm{\lambda t}}}{3})(3{\mathrm{\Sigma }}_{11}{\mu }_1{\alpha }_1+2{\mu }_1^3{\alpha }_1^2)}{{((\frac{1-e^{-2\mathrm{\lambda t}}}{2})({\mathrm{\Sigma }}_{11}+{\mu }_1^2{\alpha }_1))}^{3/2}}\end{array}$(9) ), is $\mathrm{Skew}(Y)=\frac{2^{3/2}}{3}\left(\frac{3{\sigma }^2\mu /b+2{\mu }^3/b^2}{({\sigma }^2+{\mu }^2/b{)}^{3/2}}\right).$ The results in Table 2 show how $d_{+}^{\ast }$ and $d_{-}^{\ast }$ may differ due to skewness. For negatively skewed spreads, the trade is more likely to enter from the negative level, and so the positive level is slightly lower to compensate. Figures 5 and 1 show the sample paths of the spread when $\mu =-0.5$ (large negative skew) and $\mu =0$ (no skew), respectively.

Figure 5. A sample path of the spread $X\sim \mathrm{OU}\text{-}\mathrm{VG}(1,5,-0.5,0.015,0.5)$ using the optimal entry levels in the first row of Table 2.

Table 2. For various value of μ, the optimal entry levels $d_{+}^{\ast }$ and $d_{-}^{\ast }$ are shown.

μ	$\mathrm{Skew}(Y)$	$d_{+}^{\ast }$	$d_{-}^{\ast }$
−0.5	−0.825	0.362	0.488
−0.2	−0.660	0.231	0.293
−0.05	−0.225	0.213	0.227
0	0	0.220	0.220

Notes: The other parameters are $(\lambda ,b,{\sigma }^2)=(1,5,0.015)$, $X_0=0$, $\gamma =0$, r = 0.01.

In the classical OU spread model, which has no skew, Leung and Li (2015) derive the optimal entry and exit levels with transaction costs, where the two levels can also be asymmetric.

5.3. Exit Level Different From the Mean

Recall that we set $c=c_{+}=c_{-}$ for simplicity. So far, we have set c = 0. However, it is possible to find examples where $c^{\ast }$ is much greater than 0. Such a situation can occur when $X_0>\overline{\mu }$ and the discount rate r is very large. Let $X_0=0.25$ and r = 1 while the other parameters are as before with $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,5,0,0.015,0)$. We optimize c, d for 0<c<d, and find that the optimal exit level is $c^{\ast }=0.105$, while the optimal entry level $d^{\ast }$ is any value $c^{\ast }<d^{\ast }<0.25$ as shown in Figure 6, which plots the estimate of the value function $V(c,d)$.

Figure 6. The estimate of the value function $V(c,d)$ with entry level d and exit level c. The parameters are $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,5,0,0.015,0)$, $X_0=0.25$, $\gamma =0$, r = 1.

Thus, the pairs trading strategy is to enter the trade immediately, and because the large discount rate greatly reduces the profit if the trade is not exited quickly, the trade is exited at $c^{\ast }=0.105$ above the stationary mean $\overline{\mu }$ rather than waiting for it to fall to $\overline{\mu }$.

Figure 7 shows how the optimal exit level $c^{\ast }$ changes for different values of the discount rate r. For the plotted values, the point at which $c^{\ast }>0.000$ (to 3 decimal places) is r = 0.3. For all these plotted values of r, if the initial value is instead set to $X_0=\overline{\mu }=0$, then the optimal exit level would be $c^{\ast }=0.000$. This demonstrates how having an initial value different from the stationary mean and a large discount rate can make the optimal exit level $c^{\ast }$ different from 0. It is difficult to find other situations where this occurs, which suggests optimizing c may have little benefit beyond this, and fixing c = 0, which corresponds to exiting the trade when the spread passes $\overline{\mu }$, simplifies the optimization.

Figure 7. The optimal exit level $c^{\ast }$ for various value of discount rate $r=0.01,0.1,0.2,\dots ,1.5$. The other parameters are $(\lambda ,b,\mu ,{\sigma }^2,\eta )=(1,5,0,0.015,0)$, $X_0=0.25$, $\gamma =0$.

5.4. Effect of Jumps

Now we consider the effect of the parameter b, which affects the jumps of the LDOUP X, on the optimal entry level $d^{\ast }$. As before, let $(\lambda ,\mu ,{\sigma }^2,\eta )=(1,0,0.015,0)$. As $b\to \mathrm{\infty }$, the OU-VG process for the spread converges in law to a classical OU process driven by Brownian motion, and the sample paths are closer to continuous. As $b\to 0$, roughly speaking, the sample paths become more jumpy (specifically, for each t, the BDLP $Z(t)$ converges in distribution to a constant while its kurtosis diverges). Define the overshoot for the entry level d as the random variable $O=\{\begin{array}{ll}|X({\tau }_d)-(\overline{\mu }+d)|& \mathrm{if}{\tau }_d<{\tau }_{-d},\\ |X({\tau }_{-d})-(\overline{\mu }-d)|& \mathrm{if}{\tau }_d\ge {\tau }_{-d}.\end{array}$ Table 3 shows the effect of the sample paths of the spread, roughly speaking, being very jumpy (b = 1), moderately jumpy (b = 5), and approximately continuous $(b=100)$. A plot of a sample path in each of these cases is given in Figures 8(a), 1 and 8(b), respectively. The results show that the more jumpy the sample path, the higher the level at which the trade is entered to compensate for overshooting. When b = 1, the mean overshoot is relatively large (0.0640) compared to the optimal entry level $d^{\ast }$ (0.246). This demonstrates the importance of accounting for jumps by using LDOUPs compared to the classical OU process. Both Figure 8 and Table 3 show that as the sample path becomes closer to continuous, the mean of the overshoots becomes smaller and the optimal expected profit ${\hat{V}}_M(d^{\ast })$ becomes close to $d^{\ast }$ as expected, although they do not converge since r is nonzero. In this example, ${\hat{V}}_M(d^{\ast })$ increases as b decreases. Note that the stationary variance remains constant, so the increased profitability is not due to changes in the variance.

Figure 8. A sample path of the spread $X\sim \mathrm{OU}\text{-}\mathrm{VG}(1,b,0,0.015,0)$, where (a) b = 1 and (b) b = 100, using the optimal entry levels in the first and third row of Table 3.

Table 3. For various values of b, the optimal entry level d, the estimate ${\hat{V}}_M(d^{\ast })$ of the optimal expected profit, and overshoot statistics are shown.

b	$d^{\ast }$	${\hat{V}}_M(d^{\ast })$	Sample mean O	Sample sd O
1	0.246	0.286	0.0640	0.0682
5	0.220	0.227	0.0282	0.0306
100	0.211	0.196	0.0086	0.0086

Notes: The other parameters are $(\lambda ,\mu ,{\sigma }^2,\eta )=(1,0,0.015,0)$, $X_0=0$, $\gamma =0$, r = 0.01.

6. Bivariate Pairs Trading

We now consider pairs trading in a bivariate setting, where the spread $\mathbf{X}=(X_1,X_2)$ is a bivariate LDOUP with stationary mean $\overline{\mathit{\mu }}=({\overline{\mu }}_1,{\overline{\mu }}_2)$. While we focus on the 2-dimensional case, the concepts here can in principle be extended to the n-dimensional case.

6.1. Paris Trading Problem

While there are many ways to set up a pairs trading strategy on multiple spreads, we use the following setup to generalize the notion of one trade cycle used in the univariate setting. A trade is entered for the first time that either spread $X_k$ is above ${\overline{\mu }}_k+d_k$ or below ${\overline{\mu }}_k-d_k$ for k = 1, 2. If a trade on spread $X_k$ is entered, then exit occurs when the spread passes ${\overline{\mu }}_k+c_k$ or ${\overline{\mu }}_k-c_k$, respectively. When a trade is entered, no further trades can be initiated, even if the other spread later passes the entry level.

More generally, it is possible to have asymmetric entry and exit levels for both spreads, where $d_{+,1}$, $d_{-,1}$ are the levels for entering a trade on spread $X_1$ from the positive and negative sides, respectively, and similarly for the other entry and exit levels. For each spread k = 1, 2, there are 8 relevant passage times defined by $\begin{array}{rl}{\tau }_{d_{+,k}}& =inf\{t\ge 0:X_k(t)>{\overline{\mu }}_k+d_{+,k}\}\wedge T,\\ {\tau }_{c_{+,k}}& =inf\{t\ge {\tau }_{d_{+,k}}:X_k(t)<{\overline{\mu }}_k+c_{+,k}\}\wedge T,\\ {\tau }_{d_{-,k}}& =inf\{t\ge 0:X_k(t)<{\overline{\mu }}_k-d_{-,k}\}\wedge T,\\ {\tau }_{c_{-,k}}& =inf\{t\ge {\tau }_{d_{-,k}}:X_k(t)>{\overline{\mu }}_k-c_{-,k}\}\wedge T,\end{array}$ with $0\le c_{\pm ,k}<d_{\pm ,k}$. However, we will not need to use this much generality in the examples below.

The nonzero probability of common jumps for the WVAG process used as the BDLP, and hence for $\mathbf{X}$, means that the event $\{{\tau }_{d_{+,1}}={\tau }_{d_{+,2}}\}$ also has nonzero probability, that is, both spreads can pass their corresponding entry levels simultaneously. In this case, we invest with 50% weight in each of the spreads. An alternative approach would be to invest fully in the spread furthest from the mean-reverting level. However, in the numerical examples we consider below, both approaches produce virtually the same results, since it is rare that both spreads pass the optimal entry levels simultaneously. Therefore, we will only consider the former approach.

Let $\tau =min\{{\tau }_{d_{+,1}},{\tau }_{d_{-,1}},{\tau }_{d_{+,2}},{\tau }_{d_{-,2}}\}$, $\mathbf{c}=(c_{+,1},c_{-,1},c_{+,2},c_{-,2})$ and $\mathbf{d}=(d_{+,1},d_{-,1},d_{+,2},d_{-,2})$. The profit over one trade cycle is given by $\begin{array}{rl}P:=P(\mathbf{c},\mathbf{d})& ={\mathbf{1}}_{\{\tau \in \{d_{+,1},d_{+,2}\}\}}\sum _{k=1,2}{e}^{-r{\tau }_{c_{+,k}}}{W}_k(X_k({\tau }_{d_{+,k}})-X_k({\tau }_{c_{+,k}}))\\ & +{\mathbf{1}}_{\{\tau \in \{d_{-,1},d_{-,2}\}\}}\sum _{k=1,2}{e}^{-r{\tau }_{c_{-,k}}}{W}_k(X_k({\tau }_{c_{-,k}})-X_k({\tau }_{d_{-,k}})),\end{array}$ where the weights are $(W_1,W_2)=\{\begin{array}{ll}(1,0)& \mathrm{if}{\tau }_{d_{+,1}}<{\tau }_{d_{+,2}}\mathrm{or}{\tau }_{d_{-,1}}<{\tau }_{d_{-,2}},\\ (0,1)& \mathrm{if}{\tau }_{d_{+,1}}>{\tau }_{d_{+,2}}\mathrm{or}{\tau }_{d_{-,1}}>{\tau }_{d_{-,2}},\\ (0.5,0.5)& \mathrm{otherwise},\end{array}$ and r is the discount rate.

The objective function is the value function $V(\mathbf{c},\mathbf{d})=\mathbb{E}[P]-\mathrm{\gamma Var}(P),\gamma \ge 0.$ We want to find the optimal level $\mathbf{c},\mathbf{d}$ for the pairs trading strategy to maximize the value function, that is $\underset{\mathbf{c},\mathbf{d}}{max}V(\mathbf{c},\mathbf{d}).$

6.2. Simulating the Spread

Now we specifically let the spread be $\mathbf{X}\sim \mathrm{OU}\text{-}\mathrm{WVA}{G}^2(a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$. To simulate $\mathbf{X}(t_0),\dots ,\mathbf{X}(t_q)$, we use the Euler scheme approximation in Lu (2022, Sections 4.2 and 5.1), which is based on taking a stochastic integral representation of ${\mathbf{Z}}^{\ast }(\mathrm{\Delta })$, which can be approximately simulated whenever $\mathbf{Z}(t)$ can be simulated, in this case using (5(5) $\mathbf{Z}\stackrel{D}{=}\mathit{\eta I}+{\mathbf{V}}_0+(V_1,\dots ,V_n),$(5) ). This is in contrast to the univariate OU-VG process which was simulated exactly in Section 4 with a method that has no multivariate generalization. We now determine the optimal trading levels using Monte Carlo methods.

6.3. Numerical Implementation and Analysis

We consider two examples of bivariate pairs trading on the spread $\mathbf{X}\sim \mathrm{OU}\text{-}\mathrm{WVA}{G}^2(\lambda ,a,\mathit{\alpha },\mathit{\mu },\mathrm{\Sigma },\mathit{\eta })$ simulated using the approximate method above with stepsizes $\mathrm{\Delta }=0.01$ and $\tilde{\mathrm{\Delta }}=0.001$, and terminal time T = 50. In all cases, the initial value is ${\mathbf{X}}_0=\mathbf{0}$, $\gamma =0$, we do not consider control variates, and $\mathit{\mu }=-\mathit{\eta }$ so that $\overline{\mathit{\mu }}=\mathbf{0}$. Unless stated otherwise, we assume for simplicity that the entry levels are symmetric with $d_k=d_{+,k}=d_{-,k}$, k = 1, 2, and the exit levels are $\mathbf{c}=\mathbf{0}$, so we optimize the value function over $d_1,d_2$. We use m = 10,000 Monte Carlo simulations.

In our numerical examples, we explore the effect of correlation on the optimal entry level and the optimal expected profit. In the first example, the correlation has little effect, whereas in the second example, it has a larger effect. The correlation of the components of $\mathbf{X}$ is in part controlled by the correlation parameter $\rho ={\mathrm{\Sigma }}_{12}/\sqrt{{\mathrm{\Sigma }}_{11}{\mathrm{\Sigma }}_{22}}$. However, note that ρ is the correlation of the Brownian motion subordinate in (6(6) $\mathbf{Z}\stackrel{D}{=}\mathbf{B}\odot \mathbf{T},$(6) ), and is related to the correlation of the components of $\mathbf{X}$ by (15) $\begin{array}{rl}\mathrm{Corr}(X_1(t),X_2(t))& =\mathrm{Corr}(Z_1(t),Z_2(t))\\ & =\frac{a(({\alpha }_1\wedge {\alpha }_2){\mathrm{\Sigma }}_{12}+{\alpha }_1{\alpha }_2{\mu }_1{\mu }_2)}{\sqrt{{\mathrm{\Sigma }}_{11}+{\alpha }_1{\mu }_1^2}\sqrt{{\mathrm{\Sigma }}_{22}+{\alpha }_2{\mu }_2^2}}\end{array}$(15) for all t>0, due to (10(10) $\begin{array}{rl}\mathrm{Cov}(X_1(t),X_2(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)a(({\alpha }_1\wedge {\alpha }_2){\mathrm{\Sigma }}_{12}+{\alpha }_1{\alpha }_2{\mu }_1{\mu }_2).\end{array}$(10) ) and (8(8) $\begin{array}{rl}\mathrm{Var}(X_1(t))& =\left(\frac{1-e^{-2\mathrm{\lambda t}}}{2}\right)({\mathrm{\Sigma }}_{11}+{\alpha }_1{\mu }_1^2),\end{array}$(8) ).

Example 6.1

Let (16) $\begin{array}{rl}\lambda & =1,a=2.5,\mathit{\alpha }=\left(\begin{array}{c}0.2\\ 0.3\end{array}\right),\mathit{\mu }=\left(\begin{array}{c}0\\ -0.2\end{array}\right),\mathrm{\Sigma }=\left(\begin{array}{cc}0.015& \rho \sqrt{{\mathrm{\Sigma }}_{11}{\mathrm{\Sigma }}_{22}}\\ \ast & 0.02\end{array}\right),\end{array}$(16) and r = 0.01.

Table 4 shows the optimal entry levels $d_1^{\ast }$ and $d_2^{\ast }$, the optimal expected profit ${\hat{V}}_M(d_1^{\ast },d_2^{\ast })$, and the estimated probability of which spread is traded, and each of these quantities is approximately equal across the values of ρ considered.

Table 4. For various values of ρ, the optimal entry levels $d_1^{\ast },d_2^{\ast }$, the estimate ${\hat{V}}_M(d_1^{\ast },d_2^{\ast })$ of the optimal expected profit, and the estimated probability of which spread $X_1,X_2$ is traded are shown.

ρ	0.9	0.3	0	−0.3	−0.9
$\mathrm{Corr}(X_1(t),X_2(t))$	0.356	0.119	0	−0.119	−0.356
${\hat{V}}_M(d_1^{\ast },d_2^{\ast })$	0.334	0.335	0.336	0.336	0.334
$d_1^{\ast }$	0.318	0.311	0.305	0.316	0.316
$d_2^{\ast }$	0.338	0.345	0.343	0.334	0.345
Trades only $X_1$ (%)	11.42	14.15	15.88	12.72	12.40
Trades only $X_2$ (%)	83.08	80.34	79.46	83.15	81.38
Trades both (%)	0.97	0.76	0.58	0.56	0.74
Trades neither (%)	4.53	4.75	4.08	3.57	5.48

Notes: The other parameters are given in (16(16) $\begin{array}{rl}\lambda & =1,a=2.5,\mathit{\alpha }=\left(\begin{array}{c}0.2\\ 0.3\end{array}\right),\mathit{\mu }=\left(\begin{array}{c}0\\ -0.2\end{array}\right),\mathrm{\Sigma }=\left(\begin{array}{cc}0.015& \rho \sqrt{{\mathrm{\Sigma }}_{11}{\mathrm{\Sigma }}_{22}}\\ \ast & 0.02\end{array}\right),\end{array}$(16) ) with ${\mathbf{X}}_0=\mathbf{0}$, $\gamma =0$, r = 0.01.

When pairs trading on the univariate spreads $X_1$ and $X_2$, the optimal expected profit is 0.227 (see the second row of Table 3) and 0.331 respectively, while it is 0.334–0.336 in this bivariate pairs trading example. Bivariate pairs trading is more profitable than univariate pairs trading on either spread, although it is approximately as profitable as trading only $X_2$ since this accounts for 79–84% of trades. The optimal expected profit of the bivariate pairs trade must be at least that of the best univariate pairs trade, as the latter strategy can be implemented within bivariate pairs trading by setting the other spread to have an arbitrarily large optimal entry level. Indeed, this is happening to some extent here as the optimal entry level of the spread $X_1$ has increased from 0.220 in the univariate case to 0.305–0.318 in the bivariate case. Thus, this decreases the probability of only $X_1$ entering a trade from 96.79% to 11–16% in Table 4 (the latter accounts for the probability that a trade on $X_1$ fails to be entered because a trade on $X_2$ is entered first, while the former is univariate pairs trading and does not include this). So in this example, the method has essentially identified that $X_2$ is the more profitable spread, regardless of correlation, and has accordingly increased the optimal entry level of $X_1$ to reduce the instances where it is traded.

Example 6.2

Now we consider an example where different values of the correlation parameter ρ can have a larger effect on the optimal entry level and optimal expected profit. Suppose the parameters are now (17) $\lambda =1,a=6.65,\mathit{\alpha }=\left(\begin{array}{c}0.15\\ 0.15\end{array}\right),\mathit{\mu }=\mathbf{0},\mathrm{\Sigma }=\left(\begin{array}{cc}0.015& \rho \sqrt{{\mathrm{\Sigma }}_{11}{\mathrm{\Sigma }}_{22}}\\ \ast & 0.015\end{array}\right),$(17) and r = 1.

Recall the parameter boundary $a<\frac{1}{{\alpha }_1}\wedge \frac{1}{{\alpha }_2}$. In Example 6.1, the upper bound of $\mathrm{Corr}(X_1(t),X_2(t))$ as a function of $\rho \in (-1,1)$ is 0.395 by (15(15) $\begin{array}{rl}\mathrm{Corr}(X_1(t),X_2(t))& =\mathrm{Corr}(Z_1(t),Z_2(t))\\ & =\frac{a(({\alpha }_1\wedge {\alpha }_2){\mathrm{\Sigma }}_{12}+{\alpha }_1{\alpha }_2{\mu }_1{\mu }_2)}{\sqrt{{\mathrm{\Sigma }}_{11}+{\alpha }_1{\mu }_1^2}\sqrt{{\mathrm{\Sigma }}_{22}+{\alpha }_2{\mu }_2^2}}\end{array}$(15) ). The upper bound converges to 1 as $a\to \frac{1}{{\alpha }_k}$ while keeping ${\alpha }_k,{\mu }_k,{\mathrm{\Sigma }}_{\mathrm{kk}},{\eta }_k$ constant for all k. Since $a\approx \frac{1}{{\alpha }_1}=\frac{1}{{\alpha }_2}$ here, this example represents monitoring two spreads $X_1\stackrel{D}{=}{X}_2$ that are equal in law for each fixed ρ, with $X_1$ and $X_2$ converging pathwisely as $\rho \to 1$, since in the limit, $\mathrm{Corr}(X_1(t),X_2(t))=1$ for all t with ${\mathbf{X}}_0=\mathbf{0}$. Also, appealing to the symmetry of $\mathbf{X}(t)$ for each fixed t, all components of the optimal entry levels ${\mathbf{d}}^{\ast }$ are equal to $d^{\ast }$, say. Thus, we optimize over d, the symmetric entry level for all components.

Table 5 shows the results for this example over different values of ρ. Here, the optimal expected profit ${\hat{V}}_M(d^{\ast })$ increases from 0.032 to 0.040 as ρ decreases from 0.99 to 0, this represents a 26% increase in the optimal expected profit. The case where $\rho =0.99$ is approximately the case of univariate pairs trading, which occurs in the limit as $\rho \to 1$. Since ${\hat{V}}_M(d^{\ast })$ is larger when ρ is much less than 1, the results demonstrate a situation where bivariate pairs trading is much more profitable than what is effectively pairs trading on a single spread, and where correlation has a larger effect.

Table 5. For various value of ρ, the optimal entry level $d^{\ast }$ and the estimate ${\hat{V}}_M(d^{\ast })$ of the optimal expected profit.

ρ	0.99	0.9	0.3	0	$-0.3$	$-0.9$
$\mathrm{Corr}(X_1(t),X_2(t))$	0.988	0.898	0.299	0	$-0.299$	$-0.898$
${\hat{V}}_M(d^{\ast })$	0.032	0.035	0.039	0.040	0.040	0.034
$d^{\ast }$	0.037	0.044	0.045	0.045	0.044	0.045

Notes: The other parameters are given in (17(17) $\lambda =1,a=6.65,\mathit{\alpha }=\left(\begin{array}{c}0.15\\ 0.15\end{array}\right),\mathit{\mu }=\mathbf{0},\mathrm{\Sigma }=\left(\begin{array}{cc}0.015& \rho \sqrt{{\mathrm{\Sigma }}_{11}{\mathrm{\Sigma }}_{22}}\\ \ast & 0.015\end{array}\right),$(17) ) with ${\mathbf{X}}_0=\mathbf{0}$, $\gamma =0$, r = 1.

These results are consistent with intuition. Monitoring multiple spreads of the same law presents additional opportunities to enter into the trade sooner, which increases the optimal expected profit relative to monitoring only one spread, while having a larger discount rate r substantially rewards entering the trade sooner. Thus, in this example, using bivariate pairs trading by choosing pairs that are uncorrelated (as represented by the $\rho =0$ case) is much more profitable than either univariate pairs trading or choosing highly correlated spreads (both represented by the $\rho \to 1$ case).

7. Conclusions

We have presented a Monte Carlo framework for evaluating pairs trading strategies where the mean-reverting spread follows an LDOUP. We numerically demonstrated how different parameters affect the optimal trading level and value function and capture various aspects of the trading strategy. The framework is flexible to accommodate a wide class of LDOUPs, provided that they can be fitted and simulated, as well as different ways of trading the spread. For instance, methods for fitting and simulating univariate LDOUPs with tempered stable and normal inverse Gaussian stationary distributions are provided in Valdivieso, Schoutens, and Tuerlinckx (2009). Our proposed control variates for variance reduction, can also be applied to similar problems.

For future research, one useful direction is to incorporate transaction costs into the trading problems. To that end, Do and Faff (2012) examine the profitability of pairs trading accounting for transaction costs, and Leung and Li (2015) analyse an optimal stopping approach to pairs trading with transaction costs and stop loss, and this work could be extended to spreads following LDOUPs. Another possibility is to consider large jumps in stock prices that trigger immediate liquidation due to the need to return the stocks. This can also be incorporated into the stochastic model and value function. Additional extensions of our framework include portfolio optimization problems arising from trading multiple spreads and spread selection based on statistical and other characteristics. Theoretical results on pairs trading under LDOUP models are very limited, and additional work in this direction would be useful. However with LDOUPs, these problems are complicated and unlikely to admit analytic solutions, so simulation-based approaches can be both mathematically interesting and practically useful.

Acknowledgements

Kevin Lu thanks the Department of Applied Mathematics at the University of Washington, where the substantial majority of his work on this article was completed. We thank an anonymous referee for the useful comments which have helped improve the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix. Equivalence of control variates and linear regression

Suppose we have multiple control variate $\mathbf{C}=(C_1,\dots ,C_p)$ to estimate $\mathbb{E}[P]$, where $(\mathbf{C},P)$ is jointly simulated as we want the control variates to be correlated with P. As explained in Glasserman (2003, Section 4.1.2), the control variate estimator is ${\hat{P}}_C=\overline{Y}-{\hat{\mathit{\beta }}}_C^{\mathrm{\prime }}(\overline{\mathbf{X}}-{\mathit{\mu }}_C),$ where $\mathbf{Y}$ is the vector of m simulated values of P, $\overline{Y}$ is its sample mean, ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_p$ are the vectors of m simulated values of $C_1,\dots ,C_p$, respectively, $\overline{\mathbf{X}}$ is the vector of the sample means of ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_p$, ${\mathit{\mu }}_C=\mathbb{E}[\mathbf{C}]$, ${\hat{\mathit{\beta }}}_C=S_{\mathrm{XX}}^{-1}{S}_{\mathrm{XY}}$, $S_{\mathrm{XX}}\in {\mathbb{R}}^{p\times p}$ with the sample covariance $\widehat{\mathrm{Cov}}({\mathbf{X}}_i,{\mathbf{X}}_j)$ as its $(i,j)$-entry, and $S_{\mathrm{XY}}\in {\mathbb{R}}^p$ with the sample covariance $\widehat{\mathrm{Cov}}({\mathbf{X}}_i,\mathbf{Y})$ as its ith entry. This generalizes the method of Section 4.2.

Consider the linear regression model $\mathbf{Y}=X\mathit{\beta }+\mathit{ϵ}$, where $X\in {\mathbb{R}}^{(p+1)\times m}$ is the design matrix, which includes an intercept and the control variates $C_1,\dots ,C_p$, and $\mathit{\beta }\in {\mathbb{R}}^{p+1}$ with LSE $\hat{\mathit{\beta }}$. The prediction of this linear regression at the point $\mathbf{x}=(1,{\mathit{\mu }}_C)$ is $\hat{Y}={\hat{\mathit{\beta }}}^{\mathrm{\prime }}\mathbf{x}=\overline{Y}+{\hat{\mathit{\beta }}}_C^{\mathrm{\prime }}({\mathit{\mu }}_C-\overline{\mathbf{X}})={\hat{P}}_C,$ where the second equality follows from the centred form of linear regression by Rencher and Schaalje (2008, pg 156, Equations (7.36), (7.46)). This reference also provides that ${\hat{\mathit{\beta }}}_1={\hat{\mathit{\beta }}}_C$, where $\hat{\mathit{\beta }}=({\hat{\beta }}_0,{\hat{\mathit{\beta }}}_1)$.

Thus, the control variate estimator ${\hat{P}}_C$ is the prediction of this linear regression. This fact is well-known and mentioned in Glasserman (2003, Section 4.1.1) for a single control variate, and it is unsurprising that it is also true for multiple control variates. Nelson (1990) gives a different but equivalent statement where the regressors are $\mathbf{C}-{\mathit{\mu }}_C$.