Valuation of European Options Under an Uncertain Market Price of Volatility Risk

Abstract

We propose a model to quantify the effect of parameter uncertainty on the option price in the Heston model. More precisely, we present a Hamilton–Jacobi–Bellman framework which allows us to evaluate best and worst-case scenarios under an uncertain market price of volatility risk. For the numerical approximation, the Hamilton–Jacobi–Bellman equation is reformulated to enable the solution with a finite element method. A case study with butterfly options exhibits how the dependence of Delta on the magnitude of the uncertainty is nonlinear and highly varied across the parameter regime.

Keywords:

• Uncertain market price
• volatility risk
• Hamilton–Jacobi–Bellman equation
• finite element method
• uncertainty quantification

1. Introduction

One of the main challenges in financial mathematics is to determine the fair pricing of options. Among the simplifications underlying the original Black–Scholes (BS) model is the assumption that the volatility of the stock price is constant. Analysis of the real-life data does not support this statement and so numerous attempts were made to differently model the behaviour of volatility in time. A widely used approach was proposed in Heston (1993) which models the volatility through another, correlated stochastic process. However, Heston (1993) introduces, compared to the BS model, additional parameters such as the market price of volatility risk.

Generally in option pricing, errors in parameter estimates can lead to inconsistent results even if the underlying model were to be accurate. This inspired the approach taken in Avellaneda, Levy, and Parás (1995) where instead of trying to model and predict the behaviour of parameters, they are assumed to stay within a given tolerance interval. This allows to manage the risk by considering the worst-case scenario, leading in fact to an optimal control problem involving nonlinear PDE as shown in Karoui, Peng, and Claire Quenez (1997). In Cohen and Tegnér (2019) parameter uncertainty of the Heston model is examined by means of Hamilton–Jacobi–Bellman equations, with uncertainty in the risk free rate r, long-term mean γ and reversion level κ. María et al. (2018) use stochastic collocation of the Heston model for the uncertainty qualification of the risk-free rate, long-term mean, reversion level, initial variance and correlation. We also refer to Ma and Forsyth (2016) for a related study of two-factor uncertain volatility.

The contribution of this work is to extend the Hamilton–Jacobi–Bellman framework to quantify the effect of uncertainty in the Heston model specifically in relation to the market price of volatility risk λ. This parameter is difficult to estimate in practice. For instance, in Heston (1993) a linear scaling of the market price of volatility risk is assumed and while there is an evidence of a positive correlation (see Duarte and Jones 2007) there does not appear to be a consensus on how to estimate the scaling factor (see Wu, Zhou, and Wang 2018). In fact, some authors (see for example Ikonen and Toivanen 2009; Forsyth, Vetzal, and Zvan 1998; Kunoth, Schneider, and Wiechers 2012) assume it to be equal to 0. However, other works indicate that this assumption may not hold up in various realistic settings such as those discussed for example in Doran (2007) and Bakshi and Kapadia (2003). The attempts of evaluating the market price of volatility in different financial circumstances can be found in Doran and Ronn (2008) and Wu, Zhou, and Wang (2018) while its impact on option pricing is investigated in Duarte and Jones (2007) and Mielkie and Davison (2013). Moreover, the market price of volatility risk λ is qualitatively different from those considered in Cohen and Tegnér (2019) and María et al. (2018) as λ is not part of the SDE modelling the price of the underlying security, and instead arises in the construction of a hedging portfolio. The main numerical challenge when computing the option price and Greeks associated with the uncertain Heston model is that one has to solve fully nonlinear PDEs with mixed boundary conditions. A numerical scheme capable of this is provided in Jaroszkowski and Jensen (2021) in form of a finite element method. A crucial advantage of using a finite element method in this context is that it allows us also to ensure convergence of the gradient (Jensen 2017). This is of particular importance in option pricing where partial derivatives inform the construction of a hedging portfolio.

The outline of this chapter is as follows. In Section 2 we briefly state the uncertain Heston model and then show how it can be interpreted as a backward-in-time stochastic optimal control problem. By combining the methods of stochastic volatility and uncertain parameters we obtain a second-order nonlinear PDE modelling the worst and best case scenarios when a range of values of the market price of volatility risk is considered. In Section 3 we present a transformation of the Heston equation to the form required by the numerical scheme in Jaroszkowski and Jensen (2021). Finally, in Section 4 we present a case study of a long butterfly option whose main goal is to investigate the impact of the uncertain market price of volatility risk on the option price and its derivatives.

2. The Uncertain Heston Model

In this section, we derive an extension of the Heston model to price European options in the presence of uncertain parameters.

Consider a stock with price S and a European option with expiry time $T\ge 0$ and value V. Given a Wiener process $W_1\left(t\right)$, we model the change of the stock price with (1) $\mathrm{d}S\left(t\right)=\mu S\left(t\right)\mathrm{d}t+\sigma S\left(t\right)\mathrm{d}{W}_1\left(t\right).$(1) While in the classical Black–Scholes model the volatility σ is assumed to be constant, we follow (Heston 1993) and represent it by yet another, correlated Wiener process $W_2\left(t\right)$. Then, denoting the volatility of volatility as ξ we obtain the second stochastic differential equation (2) $\mathrm{d}v\left(t\right)=\kappa (\gamma -v(t\left)\right)\mathrm{d}t+\xi \sqrt{v\left(t\right)}\mathrm{d}{W}_2\left(t\right),$(2) where the variance $v={\sigma }^2$ is a square of the volatility σ and ξ is assumed to be a known constant. We denote the correlation coefficient between $W_1$ and $W_2$ as $\rho \in (-1,1)$. Note that (2(2) $\mathrm{d}v\left(t\right)=\kappa (\gamma -v(t\left)\right)\mathrm{d}t+\xi \sqrt{v\left(t\right)}\mathrm{d}{W}_2\left(t\right),$(2) ) is a mean-reverting process with a long-term mean equal to γ and a reversion level equal to κ. We take the drift coefficient μ in (1(1) $\mathrm{d}S\left(t\right)=\mu S\left(t\right)\mathrm{d}t+\sigma S\left(t\right)\mathrm{d}{W}_1\left(t\right).$(1) ) to be constant and generally unknown.

At this point, we summarize the underlying assumptions of the model. The dividend payouts during the lifetime of the option are set to be 0. We assume it is possible to lend and borrow any amount of a risk-free asset at a known constant interest rate r. Moreover, we are allowed to trade any amount, possibly fractional, of the stock S or an option of value $V(t,S,v)$ at any time $0\le t\le T$. We also say that the market is frictionless, which means that no such transaction generates fees. Lastly, we assume that there is no arbitrage opportunity.

It is well known (Heston 1993), also Wilmott (2006, Chapter 51), that in this setting the existence of hedging portfolios ensures that the option price V admits the parabolic equation (3) $\begin{array}{r}\frac{\mathrm{\partial }V}{\mathrm{\partial }t}+\frac{1}{2}v{S}^2\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{S}^2}+Sv\xi \rho \frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }S\mathrm{\partial }v}+\frac{1}{2}v{\xi }^2\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{v}^2}+rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}+\left(\kappa \right(\gamma -v)-\xi \lambda \sqrt{v})\frac{\mathrm{\partial }V}{\mathrm{\partial }v}=0,\end{array}$(3) for some universal function $\lambda (S,\sigma ,t)$ which is common to all options. The function λ is called the market price of volatility risk because it can be interpreted as the value which market participants assign to the volatility risk. Yet, there is no agreed way to measure this function from market prices. Thus choosing the function $\lambda (S,\sigma ,t)$ is extremely challenging from the theoretical point of view and as a practical task.

Heston (1993) proposes that it can be chosen equal to $\lambda \sigma$ with $\lambda \in \mathbb{R}$ being a scaling factor, i.e., in this special case λ has the meaning of a coefficient and not the whole market price of volatility risk function.

Since there is no agreed method of estimating the market price of volatility, one may argue that any such estimate will be burdened with inaccuracies. We propose therefore a new methodology to take this lack of knowledge about the market into account: We will model the market price of volatility risk as uncertain, borrowing concepts from Avellaneda, Levy, and Parás (1995), where the effect of uncertain volatility is examined.

More concretely, we assume that λ is an unknown parameter contained in some interval $L\subset \mathbb{R}$ and consider the set $\mathbb{L}$ of all measurable mappings from $[0,T]$ to L. For all $\mathit{\lambda }\in \mathbb{L}$, we define the linear operators ${\mathcal{L}}^{\mathit{\lambda }}$ as (4) $\begin{array}{r}-\frac{1}{2}(S^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{S}^2}+2\rho \xi vS\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }S\mathrm{\partial }v}+{\xi }^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{v}^2})-rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}-\left[\kappa \right(\gamma -v)-\xi \mathit{\lambda }(t\left)\sqrt{v}\right]\frac{\mathrm{\partial }V}{\mathrm{\partial }v}+rV.\end{array}$(4) We define the Heston equation associated with the control $\mathit{\lambda }\in \mathbb{L}$ to be (5) $-{\mathrm{\partial }}_{t}V+{\mathcal{L}}^{\mathit{\lambda }}V=0.$(5) Let us momentarily assume that $\mathit{\lambda }$ is known and given. Then, in order to obtain a well-posed problem, we need to also enforce boundary and the final time conditions. Throughout this chapter, given maturity T, we will use the following conditions (6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) (6b) $\begin{array}{rl}V(0,v,t)& =\mathrm{\Lambda }\left(0\right),\end{array}$(6b) (6c) $\begin{array}{rl}\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }S}(S,v,t)& =\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }S}\left(S\right),\end{array}$(6c) (6d) $\begin{array}{rl}-rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}(S,0,t)-\kappa \gamma \frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,0,t)+rV(S,0,t)-{\mathrm{\partial }}_{t}V(S,0,t)& =0,\end{array}$(6d) (6e) $\begin{array}{rl}\underset{v\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,v,t)& =0,\end{array}$(6e) where the given function Λ is the pay-off profile of the option.

The boundary condition (6d(6d) $\begin{array}{rl}-rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}(S,0,t)-\kappa \gamma \frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,0,t)+rV(S,0,t)-{\mathrm{\partial }}_{t}V(S,0,t)& =0,\end{array}$(6d) ) for a vanishing variance can be thought of as taking limit $v\to 0$ in the Black–Scholes equation and it is adapted directly from Heston (1993). Notice on the other hand how compared to Heston (1993) the Dirichlet condition for a large volatility was replaced by the Neumann condition (6e(6e) $\begin{array}{rl}\underset{v\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,v,t)& =0,\end{array}$(6e) ), which is more favourably from the numerical point of view. We motivate these Neumann conditions with the observation that as the volatility approaches extremely large values, the influence of its oscillations on the option price is negligible. Thus we impose that the rate of change of V in v-direction at this asymptotic boundary to be 0. In the literature, such approach was adopted for the pricing of American options in Clarke and Parrott (1999) and Ikonen and Toivanen (2008), see also Zhu and Chen (2011).

Remark 2.1

Recall that for a call option with the strike price K one would choose Λ as $\mathrm{\Lambda }\left(S\right)=max(0,S-K)$. In order to calculate the value of a long butterfly position of width $2a$ and the strike price K the choice would be $\mathrm{\Lambda }\left(S\right)=max(0,S-(K-a\left)\right)-2max(0,S-K)+max(0,S-(K+a\left)\right).$ Similarly, to consider the value of a long straddle with the strike price K one would require $\mathrm{\Lambda }\left(S\right)=max(0,S-K))-max(0,S-K)$.

In practice, the implementation of the numerical scheme will require us to truncate the domain. Generally, we choose a rectangular domain $\mathrm{\Omega }=[S_{min},S_{max}]\times [0,v_{max}]$ as in Figure 1 with $\begin{array}{rl}\mathrm{\partial }{\mathrm{\Omega }}_{R_t}& =[S_{min},S_{max}]\times \left\{0\right\},\\ \mathrm{\partial }{\mathrm{\Omega }}_{R_1}& =[S_{min},S_{max}]\times \left\{v_{max}\right\},\\ \mathrm{\partial }{\mathrm{\Omega }}_{R_2}& =\left\{S_{max}\right\}\times [0,v_{max}],\\ \mathrm{\partial }{\mathrm{\Omega }}_D& =\left\{S_{min}\right\}\times [0,v_{max}].\end{array}$ For the sake of brevity, we will also define ${\mathrm{\Phi }}^{\mathit{\lambda }}$ which we call the Heston operator: (7) $\begin{array}{rl}& {\mathrm{\Phi }}^{\mathit{\lambda }}V(t,S,v)=\\ & \{\begin{array}{ll}{\displaystyle -{\mathrm{\partial }}_{t}V(t,S,v)+{\mathcal{L}}^{\mathit{\lambda }\left(t\right)}V(t,S,v)}& \mathrm{i}\mathrm{f}(t,S,v)\in [0,T)\times \mathrm{\Omega },\\ {\displaystyle V(t,S,v)-\mathrm{\Lambda }\left(0\right)}& \mathrm{i}\mathrm{f}(t,S,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_D,\\ {\displaystyle (0,1)\cdot \mathrm{\nabla }V(t,S,v)}& \mathrm{i}\mathrm{f}(t,S,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_1},\\ {\displaystyle (1,0)\cdot \mathrm{\nabla }V(t,S,v)-\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }S}\left(S\right)}& \mathrm{i}\mathrm{f}(t,S,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2},\\ {\displaystyle -{\mathrm{\partial }}_{t}V(t,S,v)-(rS,\kappa \gamma )\cdot \mathrm{\nabla }V(t,S,v)+rV(t,S,v)}& \mathrm{i}\mathrm{f}(t,S,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_t},\\ {\displaystyle V(t,S,v)-\mathrm{\Lambda }\left(S\right)}& \mathrm{i}\mathrm{f}(t,S,v)\in \left\{T\right\}\times \overline{\mathrm{\Omega }}.\end{array}\end{array}$(7)

Figure 1. Truncated domain $\mathrm{\Omega }=[S_{min},S_{max}]\times [0,v_{max}]$.

2.1. Extremal Behaviour Under Uncertainty

Of particular interest are the highest and lowest option price which can occur with market prices of volatility risk $\mathit{\lambda }\in \mathbb{L}$. For $(t,S,v)\in [0,T]\times \overline{\mathrm{\Omega }}$ we set (8a) $\begin{array}{rl}{V}_{sup}(t,S,v)& :=\underset{\mathit{\lambda }\in \mathbb{L}}{inf}\left\{V^{\mathit{\lambda }}\right(t,S,v):{\mathrm{\Phi }}^{\mathit{\lambda }}{V}^{\mathit{\lambda }}=0\mathrm{o}\mathrm{n}[0,T]\times \overline{\mathrm{\Omega }}\},\end{array}$(8a) (8b) $\begin{array}{rl}{V}_{inf}(t,S,v)& :=\underset{\mathit{\lambda }\in \mathbb{L}}{sup}\left\{V^{\mathit{\lambda }}\right(t,S,v):{\mathrm{\Phi }}^{\mathit{\lambda }}{V}^{\mathit{\lambda }}=0\mathrm{o}\mathrm{n}[0,T]\times \overline{\mathrm{\Omega }}\}.\end{array}$(8b) We seek an alternative characterization of the functions $V_{sup}$ and $V_{inf}$ as solution of a PDE, more precisely of a Hamilton–Jacobi–Bellman (HJB) equation in combination with final time and boundary conditions. We shall focus here on the argument for $V_{sup}$; the corresponding analysis for $V_{inf}$ follows analogously. Notice also that the $sup$ in the notation of $V_{sup}$ refers to the $sup$ in the below HJB Equation (10(10) $-{\mathrm{\partial }}_t{V}_{sup}(t,S,V)+\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)=0.$(10) ), not to the $sup$ on the right-hand side of (8).

Suppose there is a minimizer of the right-hand side of (8), which we denote $\stackrel{\mathbf{ˆ}}{\mathit{\lambda }}$. Then, for $(t,S,v)\in [0,T)\times \mathrm{\Omega })$ and $V_{sup}$ sufficiently smooth, (9a) $\begin{array}{rl}{V}_{sup}(t,S,v)& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}{\mathrm{\partial }}_t{V}_{sup}(\tau ,S,v\left)\mathrm{d}\tau \right]\end{array}$(9a) (9b) $\begin{array}{rl}& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}{\mathcal{L}}^{\stackrel{\mathbf{ˆ}}{\mathit{\lambda }}\left(t\right)}{V}_{sup}(\tau ,S,v\left)\mathrm{d}\tau \right]\end{array}$(9b) (9c) $\begin{array}{rl}& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}{\mathcal{L}}^{\stackrel{\mathbf{ˆ}}{\mathit{\lambda }}\left(t\right)}{V}_{sup}(\tau ,S,v\left)\mathrm{d}\tau \right]\end{array}$(9c) (9d) $\begin{array}{rl}& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}{\mathcal{L}}^{\stackrel{\mathbf{ˆ}}{\mathit{\lambda }}\left(t\right)}{V}_{sup}(t+h,S,v\left)\mathrm{d}\tau \right]\end{array}$(9d) (9e) $\begin{array}{rl}& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t+h,S,v\left)\mathrm{d}\tau \right],\end{array}$(9e) where (9e(9e) $\begin{array}{rl}& =\underset{h\to 0}{lim}\left[V_{sup}\right(t+h,S,v)-{\int }_t^{t+h}\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t+h,S,v\left)\mathrm{d}\tau \right],\end{array}$(9e) ) follows from (8). After multiplication with $1/h$ we find that $\begin{array}{rl}\underset{h\to 0}{lim}\frac{1}{h}{\int }_t^{t+h}\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(\tau ,S,v)\mathrm{d}\tau & =\underset{h\to 0}{lim}\frac{V_{sup}(t+h,S,v)-V_{sup}(t,S,v)}{h}\\ & ={\mathrm{\partial }}_t{V}_{sup}(T,S,v).\end{array}$ Assuming that $V_{sup}$ is a classical solution so that $\underset{h\to 0}{lim}\frac{1}{h}{\int }_t^{t+h}\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(\tau ,S,v)\mathrm{d}\tau =\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)$ we arrive at (10) $-{\mathrm{\partial }}_t{V}_{sup}(t,S,V)+\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)=0.$(10) Because the boundary conditions (6a(6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) ) do not depend on λ, also $V_{sup}$ satisfies them. Hence, in summary, $V_{sup}$ solves Equation (10(10) $-{\mathrm{\partial }}_t{V}_{sup}(t,S,V)+\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)=0.$(10) ) subject to the boundary conditions (6a(6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) ).

3. Transformation of the Uncertain Heston Model

We propose to use the finite element method of Jaroszkowski and Jensen (2021) for the numerical approximation of $V_{sup}$. It is shown here that this method can capture the fully nonlinear structure of (10(10) $-{\mathrm{\partial }}_t{V}_{sup}(t,S,V)+\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)=0.$(10) ) as well as the mixed (and thus discontinuous) boundary conditions (6a(6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) ) and that it will converge under mesh refinement to the viscosity solution of the final time boundary value problem. Importantly, this finite element approach has been shown (Jensen 2017) to guarantee strong convergence in the gradient of the value function, even for a degenerately elliptic HJB operator, as is the case here.

This section will now perform the transformation of the elliptic operators ${\mathcal{L}}^{\lambda }$ to their isotropic form in order to be consistent with the framework of the numerical method formulated in Chapter (Jaroszkowski and Jensen 2021). Our first goal is to remove the S dependence of the coefficients. In order to do that we let $S={\mathrm{e}}^x$. Then $\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}=S\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }S},\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{x}^2}=S^2\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{S}^2}+\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x},\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }S\mathrm{\partial }v}=\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }x\mathrm{\partial }v}$ where the transformed domain for x is denoted ${\mathrm{\Omega }}^{\prime }$, see Figure 2. Substituting into (4(4) $\begin{array}{r}-\frac{1}{2}(S^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{S}^2}+2\rho \xi vS\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }S\mathrm{\partial }v}+{\xi }^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{v}^2})-rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}-\left[\kappa \right(\gamma -v)-\xi \mathit{\lambda }(t\left)\sqrt{v}\right]\frac{\mathrm{\partial }V}{\mathrm{\partial }v}+rV.\end{array}$(4) ) we get (11) $\begin{array}{rl}{\mathcal{L}}_1^{\lambda }{V}_{sup}(x,v,t)& :=-\frac{1}{2}v(\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{x}^2}+2\rho \xi \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }x\mathrm{\partial }v}+{\xi }^2\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{v}^2})\\ & -(r-\frac{1}{2}v)\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}-\left[\kappa \right(\gamma -v)-\xi \lambda \sqrt{v}]\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }v}+r{V}_{sup}.\end{array}$(11) Similarly, with $\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }S}=\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}\frac{\mathrm{\partial }x}{\mathrm{\partial }S}={\mathrm{e}}^{-x}\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x},\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }S}=\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }x}\frac{\mathrm{\partial }x}{\mathrm{\partial }S}={\mathrm{e}}^{-x}\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }x}$ the transformed Neumann boundary condition (6c(6c) $\begin{array}{rl}\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }S}(S,v,t)& =\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }S}\left(S\right),\end{array}$(6c) ) is (12) ${\frac{\mathrm{\partial }{V}_{sup}(x,v,t)}{\mathrm{\partial }x}|}_{(x,v)\in \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }}}=\underset{x\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x\left)\right)}{\mathrm{\partial }x}.$(12) The boundary conditions on $\mathrm{\partial }{\mathrm{\Omega }}_D^{\prime }$ and $\mathrm{\partial }{\mathrm{\Omega }}_{R_1}^{\prime }$ remain as in (6b(6b) $\begin{array}{rl}V(0,v,t)& =\mathrm{\Lambda }\left(0\right),\end{array}$(6b) ) and (6e(6e) $\begin{array}{rl}\underset{v\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,v,t)& =0,\end{array}$(6e) ), in the final time condition one substitutes ${\mathrm{e}}^x$ for S and the Robin condition on $\mathrm{\partial }{\mathrm{\Omega }}_{R_t}$ is obtained by substituting v = 0 into (11(11) $\begin{array}{rl}{\mathcal{L}}_1^{\lambda }{V}_{sup}(x,v,t)& :=-\frac{1}{2}v(\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{x}^2}+2\rho \xi \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }x\mathrm{\partial }v}+{\xi }^2\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{v}^2})\\ & -(r-\frac{1}{2}v)\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}-\left[\kappa \right(\gamma -v)-\xi \lambda \sqrt{v}]\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }v}+r{V}_{sup}.\end{array}$(11) ). The Heston operator ${\mathrm{\Phi }}_1$ in the new coordinates is then defined by (13) $\begin{array}{rl}& {\mathrm{\Phi }}_1^{\mathit{\lambda }}V(t,x,v)\\ & =\{\begin{array}{ll}{\displaystyle -{\mathrm{\partial }}_{t}V(t,x,v)+{\mathcal{L}}_1^{\mathit{\lambda }\left(t\right)}V(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times {\mathrm{\Omega }}^{\mathrm{\prime }},\\ {\displaystyle V(t,x,v)-\mathrm{\Lambda }\left(0\right)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_D^{\mathrm{\prime }},\\ {\displaystyle (0,1)\cdot \mathrm{\nabla }V(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_1}^{\mathrm{\prime }},\\ {\displaystyle (1,0)\cdot \mathrm{\nabla }V(t,x,v)-\underset{x\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x\left)\right)}{\mathrm{\partial }x}}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }},\\ {\displaystyle -{\mathrm{\partial }}_{t}V(t,x,v)-(r,\kappa \gamma )\cdot \mathrm{\nabla }V(t,x,v)+rV(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_t}^{\mathrm{\prime }},\\ {\displaystyle V(t,x,v)-\mathrm{\Lambda }\left(S\right(x\left)\right)}& \mathrm{i}\mathrm{f}(t,x,v)\in \left\{T\right\}\times \overline{{\mathrm{\Omega }}^{\mathrm{\prime }}}.\end{array}\end{array}$(13)

Figure 2. Domain ${\mathrm{\Omega }}^{\mathrm{\prime }}$ after the transformation $S={\mathrm{e}}^x$.

In order to remove the second-order mixed derivative from ${\mathcal{L}}_1^{\lambda }$ we consider the following change of variables: $y=x-\frac{\rho }{\xi }v,z=\frac{\sqrt{1-{\rho }^2}}{\xi }v,$ where $y\in \mathbb{R}$ and $z\ge 0$. The domain ${\mathrm{\Omega }}^{\mathrm{\prime }}$ is transformed into ${\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}$ whose shape, in general, depends on the parameters of the numerical experiment. It is depicted in Figure 2 with the numerical values of the case study in Section 4.

For $w(y,z,t):=V_{sup}\left(x\right(y,z),v(y,z),t)$ we have that $\begin{array}{rl}{\displaystyle \frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}}& {\displaystyle =\frac{\mathrm{\partial }w}{\mathrm{\partial }y},}\\ {\displaystyle \frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }v}}& {\displaystyle =-\frac{\rho }{\xi }\frac{\mathrm{\partial }w}{\mathrm{\partial }y}+\frac{\sqrt{1-{\rho }^2}}{\xi }\frac{\mathrm{\partial }w}{\mathrm{\partial }z},}\\ {\displaystyle \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{x}^2}}& {\displaystyle =\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}.}\end{array}\begin{array}{rl}{\displaystyle \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{v}^2}}& {\displaystyle =\frac{{\rho }^2}{{\xi }^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}-\frac{2\rho \sqrt{1-{\rho }^2}}{{\xi }^2}{\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }y\mathrm{\partial }z}}}\\ & +\frac{1-{\rho }^2}{{\xi }^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{z}^2},\\ {\displaystyle \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }x\mathrm{\partial }v}}& {\displaystyle =-\frac{\rho }{\xi }\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}+\frac{\sqrt{1-{\rho }^2}}{\xi }{\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }y\mathrm{\partial }z},}}\end{array}$ Combining those results with (11(11) $\begin{array}{rl}{\mathcal{L}}_1^{\lambda }{V}_{sup}(x,v,t)& :=-\frac{1}{2}v(\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{x}^2}+2\rho \xi \frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }x\mathrm{\partial }v}+{\xi }^2\frac{{\mathrm{\partial }}^2{V}_{sup}}{\mathrm{\partial }{v}^2})\\ & -(r-\frac{1}{2}v)\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}-\left[\kappa \right(\gamma -v)-\xi \lambda \sqrt{v}]\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }v}+r{V}_{sup}.\end{array}$(11) ) we obtain the canonical formulation of every ${\mathcal{L}}^{\mathit{\lambda }}$, $\mathit{\lambda }\in \mathbb{L}$, from (5(5) $-{\mathrm{\partial }}_{t}V+{\mathcal{L}}^{\mathit{\lambda }}V=0.$(5) ) as required: (14) $\begin{array}{rl}{\mathcal{L}}_2^{\mathit{\lambda }}w& :=\frac{-\xi \sqrt{1-{\rho }^2}}{2}z\mathrm{\Delta }w+(-r+\frac{\kappa \gamma \rho }{\xi }+\frac{\frac{1}{2}\xi -\kappa \rho }{\sqrt{1-{\rho }^2}}z-\mathit{\lambda }\rho \sqrt{\frac{\xi z}{\sqrt{1-{\rho }^2}}})\frac{\mathrm{\partial }w}{\mathrm{\partial }y}\\ & +(\frac{-\kappa \gamma \sqrt{1-{\rho }^2}}{\xi }+\kappa z+\mathit{\lambda }\sqrt{\xi z\sqrt{1-{\rho }^2}})\frac{\mathrm{\partial }w}{\mathrm{\partial }z}+rw.\end{array}$(14) Now we reformulate (6a(6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) )–(6e(6e) $\begin{array}{rl}\underset{v\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,v,t)& =0,\end{array}$(6e) ) accordingly. Since $\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }x}=\frac{\mathrm{\partial }w}{\mathrm{\partial }y}$, the Neumann boundary condition (6c(6c) $\begin{array}{rl}\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }S}(S,v,t)& =\underset{S\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }}{\mathrm{\partial }S}\left(S\right),\end{array}$(6c) ) for large stock prices is obtained by simply substituting y and z into (13(13) $\begin{array}{rl}& {\mathrm{\Phi }}_1^{\mathit{\lambda }}V(t,x,v)\\ & =\{\begin{array}{ll}{\displaystyle -{\mathrm{\partial }}_{t}V(t,x,v)+{\mathcal{L}}_1^{\mathit{\lambda }\left(t\right)}V(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times {\mathrm{\Omega }}^{\mathrm{\prime }},\\ {\displaystyle V(t,x,v)-\mathrm{\Lambda }\left(0\right)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_D^{\mathrm{\prime }},\\ {\displaystyle (0,1)\cdot \mathrm{\nabla }V(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_1}^{\mathrm{\prime }},\\ {\displaystyle (1,0)\cdot \mathrm{\nabla }V(t,x,v)-\underset{x\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x\left)\right)}{\mathrm{\partial }x}}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }},\\ {\displaystyle -{\mathrm{\partial }}_{t}V(t,x,v)-(r,\kappa \gamma )\cdot \mathrm{\nabla }V(t,x,v)+rV(t,x,v)}& \mathrm{i}\mathrm{f}(t,x,v)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_t}^{\mathrm{\prime }},\\ {\displaystyle V(t,x,v)-\mathrm{\Lambda }\left(S\right(x\left)\right)}& \mathrm{i}\mathrm{f}(t,x,v)\in \left\{T\right\}\times \overline{{\mathrm{\Omega }}^{\mathrm{\prime }}}.\end{array}\end{array}$(13) ) on $[0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }}$ which results in (15) ${\frac{\mathrm{\partial }w}{\mathrm{\partial }y}|}_{(y,z)\in \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }\mathrm{\prime }}}=\underset{y\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x(y,z)\left)\right)}{\mathrm{\partial }y}$(15) Under the aforementioned change of variables the Dirichlet boundary condition (6b(6b) $\begin{array}{rl}V(0,v,t)& =\mathrm{\Lambda }\left(0\right),\end{array}$(6b) ) noticing that $\underset{S\to 0}{lim}y=-\mathrm{\infty }$ converts to ${w(y,z,t)|}_{(y,z)\in \mathrm{\partial }{\mathrm{\Omega }}_D^{\mathrm{\prime }\mathrm{\prime }}}=\mathrm{\Lambda }\left(0\right),$ while for the Neumann condition (6e(6e) $\begin{array}{rl}\underset{v\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }V}{\mathrm{\partial }v}(S,v,t)& =0,\end{array}$(6e) ) we use the fact that $\frac{\mathrm{\partial }{V}_{sup}}{\mathrm{\partial }v}=-\frac{\rho }{\xi }\frac{\mathrm{\partial }w}{\mathrm{\partial }y}+\frac{\sqrt{1-{\rho }^2}}{\xi }\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$ and $\underset{v\to \mathrm{\infty }}{lim}z=\mathrm{\infty }$ to obtain $(\frac{-\rho }{\xi },\frac{\sqrt{1-{\rho }^2}}{\xi })\cdot \mathrm{\nabla }w=0.$ Analogously to the result in Heston (1993), the Robin boundary condition for $v\to 0$ is obtained simply by substituting z = 0 into (14(14) $\begin{array}{rl}{\mathcal{L}}_2^{\mathit{\lambda }}w& :=\frac{-\xi \sqrt{1-{\rho }^2}}{2}z\mathrm{\Delta }w+(-r+\frac{\kappa \gamma \rho }{\xi }+\frac{\frac{1}{2}\xi -\kappa \rho }{\sqrt{1-{\rho }^2}}z-\mathit{\lambda }\rho \sqrt{\frac{\xi z}{\sqrt{1-{\rho }^2}}})\frac{\mathrm{\partial }w}{\mathrm{\partial }y}\\ & +(\frac{-\kappa \gamma \sqrt{1-{\rho }^2}}{\xi }+\kappa z+\mathit{\lambda }\sqrt{\xi z\sqrt{1-{\rho }^2}})\frac{\mathrm{\partial }w}{\mathrm{\partial }z}+rw.\end{array}$(14) ) which gives (16) $-{\mathrm{\partial }}_{t}w+(-r+\frac{\kappa \gamma \rho }{\xi })\frac{\mathrm{\partial }w}{\mathrm{\partial }y}+\left(\frac{-\kappa \gamma \sqrt{1-{\rho }^2}}{\xi }\right)\frac{\mathrm{\partial }w}{\mathrm{\partial }z}+rw=0.$(16) We summarize the above results by introducing the transformed Heston operator ${\mathrm{\Phi }}_2$ defined as follows (17) $\begin{array}{rl}& {\mathrm{\Phi }}_2^{\mathit{\lambda }}w(t,y,z)\\ & =\{\begin{array}{ll}{\displaystyle -{\mathrm{\partial }}_{t}w(t,y,z)+{\mathcal{L}}_2^{\mathit{\lambda }\left(t\right)}w(t,y,z)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times {\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle w(t,y,z)-\mathrm{\Lambda }\left(0\right)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_D^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle (\frac{-\rho }{\xi },\frac{\sqrt{1-{\rho }^2}}{\xi })\cdot \mathrm{\nabla }w(t,y,z)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_1}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle (1,0)\cdot \mathrm{\nabla }w(t,y,z)-\underset{y\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x(y,z)\left)\right)}{\mathrm{\partial }y}}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle -{\mathrm{\partial }}_{t}w(t,y,z)-(r+\frac{\kappa \gamma \rho }{\xi },\frac{-\kappa \gamma \sqrt{1-{\rho }^2}}{\xi })}\\ \cdot \mathrm{\nabla }{V}_{sup}+rw(t,y,z)& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_t}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle w(t,y,z)-\mathrm{\Lambda }\left(S\right(x(y,z)\left)\right)}& \mathrm{i}\mathrm{f}(t,y,z)\in \left\{T\right\}\times \overline{{\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}}.\end{array}\end{array}$(17) By replacing the Heston operator ${\mathrm{\Phi }}^{\mathit{\lambda }}$ from (4(4) $\begin{array}{r}-\frac{1}{2}(S^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{S}^2}+2\rho \xi vS\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }S\mathrm{\partial }v}+{\xi }^{2}v\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{v}^2})-rS\frac{\mathrm{\partial }V}{\mathrm{\partial }S}-\left[\kappa \right(\gamma -v)-\xi \mathit{\lambda }(t\left)\sqrt{v}\right]\frac{\mathrm{\partial }V}{\mathrm{\partial }v}+rV.\end{array}$(4) ) with its transformed version from (17(17) $\begin{array}{rl}& {\mathrm{\Phi }}_2^{\mathit{\lambda }}w(t,y,z)\\ & =\{\begin{array}{ll}{\displaystyle -{\mathrm{\partial }}_{t}w(t,y,z)+{\mathcal{L}}_2^{\mathit{\lambda }\left(t\right)}w(t,y,z)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times {\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle w(t,y,z)-\mathrm{\Lambda }\left(0\right)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_D^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle (\frac{-\rho }{\xi },\frac{\sqrt{1-{\rho }^2}}{\xi })\cdot \mathrm{\nabla }w(t,y,z)}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_1}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle (1,0)\cdot \mathrm{\nabla }w(t,y,z)-\underset{y\to \mathrm{\infty }}{lim}\frac{\mathrm{\partial }\mathrm{\Lambda }\left(S\right(x(y,z)\left)\right)}{\mathrm{\partial }y}}& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_2}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle -{\mathrm{\partial }}_{t}w(t,y,z)-(r+\frac{\kappa \gamma \rho }{\xi },\frac{-\kappa \gamma \sqrt{1-{\rho }^2}}{\xi })}\\ \cdot \mathrm{\nabla }{V}_{sup}+rw(t,y,z)& \mathrm{i}\mathrm{f}(t,y,z)\in [0,T)\times \mathrm{\partial }{\mathrm{\Omega }}_{R_t}^{\mathrm{\prime }\mathrm{\prime }},\\ {\displaystyle w(t,y,z)-\mathrm{\Lambda }\left(S\right(x(y,z)\left)\right)}& \mathrm{i}\mathrm{f}(t,y,z)\in \left\{T\right\}\times \overline{{\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}}.\end{array}\end{array}$(17) ) and following the same argument as in the previous section we obtain an optimal control problem (18) $\underset{\lambda \in L}{sup}{\mathrm{\Phi }}_2^{\lambda }w(t,y,z)=0\mathrm{\forall }(t,y,z)\in [0,T]\times \overline{{\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}}$(18) analogous to (6a(6a) $\begin{array}{rl}V(S,v,T)& =\mathrm{\Lambda }\left(S\right),\end{array}$(6a) ) and (10(10) $-{\mathrm{\partial }}_t{V}_{sup}(t,S,V)+\underset{\lambda \in L}{sup}{\mathcal{L}}^{\lambda }{V}_{sup}(t,S,v)=0.$(10) ) with the structure conforming to the setting of Jaroszkowski and Jensen (2021). Note that it resembles the ‘worst-case scenario’ described in Avellaneda, Levy, and Parás (1995) but with λ instead of σ taking the role of the uncertain parameter. Having completed the transformation, this allows us to treat the market price of volatility risk as a control in an isotropic HJB problem.

4. Case Study

We now investigate the qualitative and quantitative effects of the market price of volatility risk on the price of an option. The computations use the finite element discretisation of Jaroszkowski and Jensen (2021). The code is available from the public repository (Jaroszkowski 2021) under the GNU Lesser General Public License.

As a benchmark problem, we consider the parameters of the experiment described in Doran (2004, Table 11). Following this setting we choose the final time T = 0.5, the strike price K = 50, the volatility of volatility $\xi =0.7$, the long-term volatility mean $\gamma =0.3$, the mean reversion rate $\kappa =7$ and the correlation parameter $\rho =0.5$. Since the risk-free rate r is not stated explicitly, we made a choice of r = 0.03. The domain is truncated with $v\in [0,3]$ and $S\in [1,100]$, which results in the transformed domain seen in Figure 3. In the source (Doran 2004) a simple call option is studied, which for example with $L=[-2.4,-1.6]$ has the constant optimal control $\hat{\mathit{\lambda }}\equiv -2.4$. Instead, we turn our attention here to a long butterfly position of width 40 which, as mentioned in Remark 2.1, is equivalent to choosing $\mathrm{\Lambda }\left(S\right)=max(0,S-30)-2max(0,S-50)+max(0,S-70).$

Figure 3. Domain ${\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}$ with parameter values $\rho =0.5$, $\xi =0.7$

4.1. Result 1: Value Function

We let $L=[-2.4,-1.6]$. Note how interval L is centred around the market price of volatility risk equal to $-2$ used in Doran (2004). The numerical approximation of the solution to the HJB problem is performed on the transformed domain ${\mathrm{\Omega }}^{\mathrm{\prime }\mathrm{\prime }}$ and then the resulting function is cast back to original domain $\mathrm{\Omega }$. The outcome is depicted in Figure 4. Moreover, one can see in Figure 5(a) that the numerical method in fact selects different controls as optimal in different areas of the domain. The difference between the solution of the nonlinear problem compared to the solution of the linear evolution problem associated with one of the controls can be seen in Figure 5(b). This highlights the importance of using a nonlinear model.

Figure 4. Value of a long butterfly position at t = 0 with T = 0.5, K = 50 and control set $L=[-2.4,-1.6]$.

Figure 5. Measurement of the effect of nonlinearity for a long butterfly position at $t\approx 0.39$ with T = 0.5, K = 50 and control set $L=[-2.4,-1.6]$. (a) Selected optimal control; which are throughout at the extreme points of the control set, (b) difference between solutions of a nonlinear problem and linear evolution problem with a fixed control $\lambda =-2.4$.

4.2. Result 2: λ Interval Testing

We now assess the impact of different choices of control sets L on the option value estimate. Indeed, we consider control sets of increasing diameter and measure the difference between the value function $V_{inf}$ of the worst-case scenario and the value function $V_{sup}$ of the best case scenario. The results are shown in Figure 6. The computations show the significant effect of the uncertainty in the market price of volatility risk on the option price. As indicated by Figure 6(a) the option value of worse and best case scenario can differ up to $16\mathrm{\%}$. Note that in this case, the control set contains values ranging between 0 and $-2.5$, which were found to be used in the literature. Given the evidence (see for example findings in Bakshi and Kapadia (2003)) that the market price of volatility takes negative values, the simplification of taking $\lambda =0$ may lead to erroneous estimates. On the other hand, the experiments indicate a linear correlation and in general, more negative market prices of volatility risk lead to higher option values.

Figure 6. Measurement of effect of a diameter of a control set on the value function and its derivative. Control sets are symmetrical and centred at $-1.25$, measurements were made at t = 0. (a) Comparison of $V_{sup}$ and $V_{inf}$ at $(S,v)=(2.11,2.06)$, (b) comparison of $\mathrm{\partial }{V}_{sup}/\mathrm{\partial }S$ and $\mathrm{\partial }{V}_{inf}/\mathrm{\partial }S$ at $(S,v)=(53.12,0.75)$, (c) comparison of $\mathrm{\partial }{V}_{sup}/\mathrm{\partial }S$ and $\mathrm{\partial }{V}_{inf}/\mathrm{\partial }S$ at $(S,v)=(51.76,2.84)$, (d) comparison of $\mathrm{\partial }{V}_{sup}/\mathrm{\partial }S$ and $\mathrm{\partial }{V}_{inf}/\mathrm{\partial }S$ at $(S,v)=(51.43,0.23)$.

We now direct our attention to the partial derivatives of option value V since they are used to create hedging portfolios. We investigate the effect of λ on the partial derivative of the option value with respect to S. As seen in Figure 6(b–d) the impact of the value of λ on Delta $\mathrm{\partial }V/\mathrm{\partial }S$ is strongly nonlinear in the vicinity of the strike price K. We remark at this point that numerical methods which do not guarantee gradient convergence may in general fail to capture this kind of behaviour.

Figure 6 provides a visual portrayal of the sensitivity of the price and Delta on the magnitude of the uncertainty.

4.3. Result 3: Delta Plots

In line with the results of the previous experiment, we continue to investigate the worst and the best case scenarios for the control set $L=[-2.5,0.0]$ at time t = 0. We plot differences between the Deltas $\mathrm{\partial }{V}_{sup}/\mathrm{\partial }S$ and $\mathrm{\partial }{V}_{inf}/\mathrm{\partial }S$ for all points in $\mathrm{\Omega }$ at time t = 0. The results for a call option are shown in Figure 7(a) and for a long butterfly option in Figure 7(b). Note that since $\mathrm{\partial }{V}_{inf}/\mathrm{\partial }S$ and $\mathrm{\partial }{V}_{sup}/\mathrm{\partial }S$ are both of order 1, the graphs represent a relative as well as an absolute error. We conclude that the impact of the market price of volatility risk on the delta values is significant. In the covered examples, one can expect up to $6\mathrm{\%}$ difference between the scenario where λ is neglected and the one where the HJB approach is used.

Figure 7. Comparison of plots of $\delta (V_{sup}-V_{inf})/\delta S$ at time t = 0 with control set $[-2.5,0.0]$. (a) Call option, (b) butterfly option.

Disclosure statement

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

Additional information

Funding

Bartosz Jaroszkowski acknowledges the support of the Engineering and Physical Sciences Research Council (EPSRC) grant 1816514. Max Jensen acknowledges the support of the Dr Perry James Browne Research Centre.