Toward an efficient hybrid method for pricing barrier options on assets with stochastic volatility

We consider the joint evolution of the asset price, $S_{t}$, and its stochastic variance, $V_{t}$, as follows:

$$\begin{array}[]{c}\frac{dS_{t}}{S_{t}}=rdt+\sqrt{V_{t}}\left(\rho dB_{t}+\sqrt{1-\rho^{2}}dW_{t}\right),\ \ \ S_{0}=s,\\ dV_{t}=\Phi\left(V_{t}\right)dt+\Psi\left(V_{t}\right)dB_{t},\ \ \ V_{0}=v,\end{array}$$

(6)

where $B_{t}$ and $W_{t}$ are independent Brownian motions. Given that Eqs (6) are scale invariant with respect to $S_{t}$, we can write them in terms of $X_{t}=\ln\left(S_{t}/S_{0}\right)$ and $V_{t}$:

$$\begin{array}[]{c}dX_{t}=\left(r-\frac{1}{2}V_{t}\right)dt+\sqrt{V_{t}}\left(\rho dB_{t}+\sqrt{1-\rho^{2}}dW_{t}\right),\ \ \ \ X_{0}=0,\\ dV_{t}=\Phi\left(V_{t}\right)dt+\Psi\left(V_{t}\right)dB_{t},\ \ \ V_{0}=v,\end{array}$$

(7)

Alternatively, we can study the joint evolution of the price $S_{t}$ and its volatility $\sigma_{t}$:

$$\begin{array}[]{c}\frac{dS_{t}}{S_{t}}=rdt+\sigma_{t}\left(\rho dB_{t}+\sqrt{1-\rho^{2}}dW_{t}\right),\ \ \ S_{0}=s,\\ d\sigma_{t}=\phi\left(\sigma_{t}\right)dt+\psi\left(\sigma_{t}\right)dB_{t},\ \ \ \sigma_{0}=\sigma.\end{array}$$

(8)

Given that Eqs (8) are scale invariant with respect to $S_{t}$, we can write them in terms of $X_{t}=\ln\left(S_{t}/S_{0}\right)$ and $\sigma_{t}$:

$$\begin{array}[]{c}dX_{t}=\left(r-\frac{1}{2}\sigma_{t}^{2}\right)dt+\sigma_{t}\left(\rho dB_{t}+\sqrt{1-\rho^{2}}dW_{t}\right),\ \ \ X_{0}=0,\\ d\sigma_{t}=\phi\left(\sigma_{t}\right)dt+\psi\left(\sigma_{t}\right)dB_{t},\ \ \ \sigma_{0}=\sigma,\end{array}$$

(9)

which is often more convenient. From now on, we shall concentrate of studying the dynamics of $X_{t}$.

In the general case, we can write $dB_{t}$ in the form

$$dB_{t}=\frac{dV_{t}-\Phi\left(V_{t}\right)dt}{\Psi\left(V_{t}\right)},$$

(10)

and obtain the following conditionally-independent dynamics for the log-price $X_{t}$:

$$dX_{t}=\left(r-\frac{V_{t}}{2}-\frac{\rho\sqrt{V_{t}}\Phi\left(V_{t}\right)}{\Psi\left(V_{t}\right)}\right)dt+\frac{\rho\sqrt{V_{t}}dV_{t}}{\Psi\left(V_{t}\right)}+\sqrt{1-\rho^{2}}\sqrt{V_{t}}dW_{t},\ \ \ X_{0}=0.$$

(11)

Similarly, we can write $dB_{t}$ in the form

$$dB_{t}=\frac{d\sigma_{t}-\phi\left(\sigma_{t}\right)}{\psi\left(\sigma_{t}\right)},$$

(12)

and get the following dynamics for $X_{t}$:

$$dX_{t}=\left(r-\frac{1}{2}\sigma_{t}^{2}-\frac{\rho\sigma\phi\left(\sigma_{t}\right)}{\psi\left(\sigma_{t}\right)}\right)dt+\frac{\rho\sigma_{t}d\sigma_{t}}{\psi\left(\sigma_{t}\right)}+\sqrt{1-\rho^{2}}\sigma_{t}dW_{t},\ \ \ X_{0}=0.$$

(13)

Assuming that the variance or volatility paths are given, Eqs (11), (13) describe drifted arithmetic Brownian motion with time-dependent drift and volatility.

For the well-known Heston model, Heston (1993), we have

$$\begin{array}[]{c}\Phi\left(V_{t}\right)=\kappa\left(\theta-V_{t}\right),\ \ \ \Psi\left(V_{t}\right)=\varepsilon\sqrt{V_{t}},\\ dV_{t}=\kappa\left(\theta-V_{t}\right)dt+\varepsilon\sqrt{V_{t}}dB_{t},\end{array}$$

(14)

so that Eq. (11) has the form

$$dX_{t}=\left(r-\frac{\rho\kappa\theta}{\varepsilon}-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)V_{t}\right)dt+\frac{\rho}{\varepsilon}dV_{t}+\sqrt{1-\rho^{2}}\sqrt{V_{t}}dW_{t}.$$

(15)

Thus, $X_{t}$ is the so-called drifted Brownian motion driven by the stochastic differential equation of the form

$$\begin{array}[]{c}dX_{t}=\mu\left(t\right)dt+\nu\left(t\right)dW_{t},\\ \mu\left(t\right)=\left(r-\frac{\rho\kappa\theta}{\varepsilon}-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)V_{t}\right)+\frac{\rho}{\varepsilon}\frac{dV_{t}}{dt},\\ \nu\left(t\right)=\sqrt{1-\rho^{2}}\sqrt{V_{t}}.\end{array}$$

(16)

For the Stein-Stein model, Stein and Stein (1991), we have

$$\begin{array}[]{c}\phi\left(\sigma_{t}\right)=\hat{\kappa}\left(\hat{\theta}-\sigma_{t}\right),\ \ \ \psi\left(\sigma_{t}\right)=\hat{\varepsilon},\\ d\sigma_{t}=\hat{\kappa}\left(\hat{\theta}-\sigma_{t}\right)dt+\hat{\varepsilon}dB_{t},\end{array}$$

(17)

so that Eq. (13) becomes

$$dX_{t}=\left(r-\frac{\rho\hat{\kappa}\hat{\theta}}{\hat{\varepsilon}}\sigma_{t}-\left(\frac{1}{2}-\frac{\rho\hat{\kappa}}{\hat{\varepsilon}}\right)\sigma_{t}^{2}\right)dt+\frac{\rho}{\hat{\varepsilon}}\sigma_{t}d\sigma_{t}+\sqrt{1-\rho^{2}}\sigma_{t}dW_{t}.$$

(18)

Traditionally, the Heston model is viewed as better describing the market than the Stein-Stein model since the latter does allow for zero variance. In our opinion, this is not a particularly important issue, outweighed by many advantages of the Stein-Stein model, such as a very straightforward way of simulating the evolution of the volatility path. The Heston model gained popularity due to the simple fact that it is exactly solvable for vanilla options. The explicit solution can be calculated via the original Heston (1993) formula. However, the Lewis-Lipton formula is much more efficient; see Lewis (2000), Lipton (2001), Lipton (2002), Lipton and Sepp (2008), and Schmeltze (2010).

In this paper, we consider the Heston model with constant coefficients to follow a long-established tradition, even though it is not necessarily our preferred model. We emphasize that our method is general and can handle any reasonable stochastic volatility model.

The popularity of the Heston model stems from the fact that one can write its solution in the closed-form; see Heston (1993). However, experience has shown that using the original formula is difficult due to several technical drawbacks. Therefore, for benchmarking purposes, here we present the Lewis-Lipton formula, which is easy to implement and use in practice:

$$\begin{array}[]{c}C^{H}\left(0,S_{0},K,T;r,\rho,\kappa,\theta,\varepsilon,V_{0}\right)\\ =S_{0}\left(1-\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\frac{e^{\left(i\chi+1/2\right)\left(\ln\left(K/S_{0}\right)-rT\right)+\alpha\left(T,\chi\right)-\left(\chi^{2}+1/4\right)\beta\left(T,\chi\right)V_{0}}}{\left(\chi^{2}+1/4\right)}d\chi\right),\\ \alpha\left(T,\chi\right)=-\frac{\kappa\theta}{\varepsilon^{2}}\left[\psi_{+}T+2\ln\left(\frac{\psi_{-}+\psi_{+}\exp\left(-\zeta T\right)}{2\zeta}\right)\right],\\ \beta\left(T,\chi\right)=\frac{1-\exp\left(-\zeta T\right)}{\psi_{-}+\psi_{+}\exp\left(-\zeta T\right)},\\ \psi_{\pm}=\mp\left(i\rho\varepsilon\chi+\hat{\kappa}\right)+\zeta,\\ \zeta=\sqrt{\varepsilon^{2}\left(1-\rho^{2}\right)\chi^{2}+2i\varepsilon\rho\hat{\kappa}\chi+\hat{\kappa}^{2}+\frac{\varepsilon^{2}}{4}},\end{array}$$

(19)

where $\hat{\kappa}=\kappa-\rho\varepsilon/2$. Further details are given in Lewis (2000), Lipton (2001), Lipton (2002), and Schmeltze (2010).²²2There is a typo in Lipton (2002) - a minus sign in front of $\beta$. This typo is corrected in Lipton (2018), Chapter 10. It is clear that Eq. (19) is a generalization of Eq. (4).

Unfortunately, with very few exceptions, finding a closed-form solution for barrier or other exotic options on assets with stochastic volatility is not possible, even if such a solution exists for vanilla options. Hence, more general volatility models for barrier options are as good (or bad) as the more traditional Heston and Stein-Stein models, which enjoy closed-form solutions for vanilla options.

We express the log-return process as a linear combination of the two processes:

$$X_{t}=Y_{t}+\left(\left(r-\frac{\rho\kappa\theta}{\varepsilon}\right)t-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)I_{t}\right)+\frac{\rho}{\varepsilon}\left(V_{t}-V_{0}\right)\equiv Y_{t}+M_{t},$$

(20)

where

$$\begin{array}[]{c}Y_{t}=\sqrt{1-\rho^{2}}\int_{0}^{t}\sqrt{V_{t}}dW_{t},\ Y_{0}=0,\\ I_{t}=\int_{0}^{t}V_{t^{\prime}}dt^{\prime},\ I_{0}=0,\\ M_{t}=\left(\left(r-\frac{\rho\kappa\theta}{\varepsilon}\right)t-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)I_{t}\right)+\frac{\rho}{\varepsilon}\left(V_{t}-V_{0}\right),\ \ \ M_{0}=0.\end{array}$$

(21)

Accordingly, conditionally on the filtration generated by the variance process $V_{t}$, we can represent the solution to the price process given by Eq. (6) as follows:

$$S_{t}=e^{M_{t}+Y_{t}}S_{0},$$

(22)

where $M_{t}$ is interpreted as a time-deterministic cumulative drift, and $Y_{t}$ is a martingale with a deterministic time-dependent quadratic variance.

It is clear that pricing for path-independent options, such as European calls and puts, simplifies to the Black-Scholes-Merton formula, provided that either a variance path $\left\{\left.V_{t}\right|0\leq t\leq T\right\}$ (or a volatility path $\left\{\left.\sigma_{t}\right|0\leq t\leq T\right\}$) is given. To be concrete, consider a call option with maturity $T$ and strike $K$. Eq. (6) yields

$$\frac{dS_{t}}{S_{t}}=\left(\left(r-\frac{1}{2}V_{t}\right)dt+\rho\sqrt{V_{t}}dB_{t}\right)+\sqrt{1-\rho^{2}}\sqrt{V_{t}}dW_{t},\ \ \ ,\ \ S_{0}=s.$$

(23)

Thus,

$$S_{T}=e^{rT-\frac{1}{2}\left(1-\rho^{2}\right)I_{T}+\sqrt{1-\rho^{2}}\sqrt{I_{T}}\eta}\left(e^{-\frac{1}{2}\rho^{2}I_{T}+\rho J_{T}}S_{0}\right),$$

(24)

where the non-dimensional random variables $I_{T},J_{T}$, are given by

$$I_{T}=\int_{0}^{T}V_{t}dt,\ \ \ J_{T}=\int_{0}^{T}\sqrt{V_{t}}dB_{t},$$

(25)

and $\eta$ is the standard $\left(0,1\right)$ normal variable. Accordingly for a particular trajectory $\left\{\left.V_{t}\right|0\leq t\leq T\right\}$, we obtain the following expression

$$C=C^{BS}\left(0,e^{-\frac{1}{2}\rho^{2}I_{T}+\rho J_{T}}S_{0},T,K;r,\sqrt{1-\rho^{2}}\sqrt{\frac{I_{T}}{T}}\right),$$

(26)

where the values $\left(I_{T},J_{T}\right)$ are assumed to be known. The unconditional price is obtained by averaging over all possible $\left(I_{T},J_{T}\right)$:

$$\begin{array}[]{c}C^{H}\left(0,S_{0},K,T;r,\rho,\kappa,\theta,\varepsilon,v_{0}\right)\\ =\int_{0}^{\infty}\int_{0}^{\infty}C^{BS}\left(0,e^{-\frac{1}{2}\rho^{2}I_{T}+\rho J_{T}}S_{0},T,K;r,\sqrt{1-\rho^{2}}\sqrt{\frac{I}{T}}\right)\Xi\left(I_{T},J_{T}\right)dJ_{T}dI_{T},\end{array}$$

(27)

where $\Xi\left(I_{T},J_{T}\right)$ is the joint probability density function (pdf) for the pair $\left(I_{T},J_{T}\right)$; see Willard (1997) and Romano and Touzi (1997). Thus, Eq. (27) splits the calculation of the call option price into two stages. First, the conditional price is found analytically via the standard Black-Scholes formula. Second, this conditional price is averaged according to the particular choice of the process for the variance $V_{t}$. Of course, the first stage is trivial. The usefulness of this formula depends on how easy (or difficult) it is to find the pdf for $\left(I,J\right)$ and complete the second stage.

Two approaches have been used in practice - the standard Monte-Carlo method for calculating $\left\{\left.V_{t}\right|0\leq t\leq T\right\}$, $I$, $J$, and a more advanced (but much harder) method based on the augmentation principle described in Section (13.2) Lipton (2001). Specifically, the augmented dynamic equation for $V_{t}$ yields the following system of degenerate PDEs for the triple $\left(V,I,J\right)$:

$$\begin{array}[]{c}dV_{t}=\Phi\left(V_{t}\right)dt+\Psi\left(V_{t}\right)dB_{t},\ \ \ V_{0}=v,\\ dI_{t}=V_{t}dt,\ \ \ I_{0}=0,\\ dJ_{t}=\sqrt{V_{t}}dB_{t},\ \ \ J_{0}=0.\end{array}$$

(28)

The corresponding Green’s function $G\left(t,V,I,J\right)$ is governed by the degenerate Fokker-Planck equation of the form

$$\begin{array}[]{c}G_{t}+\left(\Phi\left(V\right)G\right)_{V}+VG_{I}-\frac{1}{2}\left(\Psi^{2}\left(V\right)G\right)_{VV}-\left(\sqrt{V}\Psi\left(V\right)G\right)_{VJ}-\frac{1}{2}\left(VG\right)_{JJ}=0,\\ G\left(0,V,I,J\right)=\delta\left(V-v\right)\delta\left(I\right)\delta\left(J\right).\end{array}$$

(29)

In general, solving this equation is complicated; however, for the Heston model, Lipton (2001) found a closed-form solution.

In Figure 1 we compare prices of European call options given by analytical formula (19), and Monte Carlo formula (27). As expected, both formulas agree, although the latter formula is much slower than the former.

Our task is to price a barrier option written on an asset with stochastic volatility. For brevity, we consider barrier options with the lower barrier $B<S_{0}$ only. Considering other possibilities, such as pricing options with the upper barrier or popular double-no-touch options, is left to the reader. The corresponding IBVP can be written in the form

$$\begin{array}[]{c}P_{t}+\left(r-\frac{1}{2}V\right)P_{X}+\kappa\left(\theta-V\right)P_{V}+\frac{1}{2}V\left(P_{XX}+2\rho P_{XV}+P_{VV}\right)-rP=0,\\ P\left(T,X,V\right)=\Pi\left(X\right),\\ P\left(t,\xi,V\right)=0,\end{array}$$

(30)

where $\xi=\ln\left(B/S_{0}\right)<0$, $\Pi(X)$ is the terminal payoff function. Typical examples include the no-touch, call and put payoffs defined by:

$$\Pi\left(X_{T}\right)=1,\ \Pi\left(X_{T}\right)=S_{0}\left(e^{X_{T}}-e^{k}\right)_{+},\ \Pi\left(X_{T}\right)=S_{0}\left(e^{k}-e^{X_{T}}\right)_{+},$$

(31)

where $k=\ln\left(K/S_{0}\right)$. Alternatively, we can write the adjoint problem for Green’s function $G$:

$$\begin{array}[]{c}G_{t}+\left(\left(r-\frac{1}{2}V\right)G\right)_{X}+\left(\kappa\left(\theta-V\right)G\right)_{V}\\ -\frac{1}{2}\left(\left(VG\right)_{XX}+2\rho\left(VG\right)_{XV}+\left(VG\right)_{VV}\right)+rG=0,\\ G\left(0,X,V\right)=\delta\left(X\right)\delta\left(V-V_{0}\right),\\ G\left(t,\xi,V\right)=0,\end{array}$$

(32)

Once the corresponding Green’s function is calculated, we can find $P$ via simple integration:

$$P\left(0,0,V_{0}\right)=\int_{\xi}^{\infty}\int_{0}^{\infty}G\left(0,0,V_{0},T,X_{T},V_{T}\right)\Pi\left(X_{T}\right)dV_{T}dX_{T}.$$

(33)

While such options can be priced via either FDM or MCS, both are notoriously slow. Therefore, we want to design a much faster method, enjoying equal or higher accuracy than the classical alternatives.

As far as analytical solutions are concerned, only one is known. It was discovered by Lipton (1997), see also Lipton (2001), Lipton and McGhee (2002), and Andreasen (2001). Lipton (1997) observed that in the special case $r=0$, $\rho=0$, IBVPs (30), (32) are symmetric with respect to the transformation $X\rightarrow-X$. Hence, the classical method of images is applicable, and solutions to these problems can be presented as a linear combination of solutions without barriers. Of course, one can use this approach for options in the presence of an upper barrier, as well as for double-barrier options.

Recently, there were several unsuccessful attempts to solve the pricing problem with $r^{2}+\rho^{2}>0$. For example, De Gennaro Aquino and Bernard (2019) presented a solution, relying on an explicit expression for the joint distribution of the value of a Brownian motion with time-dependent drift and its maximum and minimum; it was quickly shown by one of the present authors that their calculation is erroneous; see De Gennaro Aquino and Bernard (2021).³³3Finding the joint distribution for a Brownian motion with time-dependent drift and its maximum and minimum is challenging. We present its solution in Section 5. He and Lin (2021) presented a “solution”, which relies on the unsubstantiated replacement of the time-dependent drift by a constant. Their approach is so arbitrary and frivolous that its detailed repudiation is not warranted.

It is hard to extend interesting formula (27) for barrier options. However, it is not impossible! Following Section13.3 in Lipton (2001), we express the value of a path-dependent option as an integral in the functional space of price trajectories:

$$P(S_{0})=e^{-rT}\int_{\Omega}\mathcal{F}(\omega)d\mathcal{D}(\omega)$$

(34)

where $\mathcal{F}(\omega)$ is a functional mapping of the space of trajectories into payoffs, and $\mathcal{D}(\omega)$ is the risk-neutral probability measure.

Further, by applying the augmentation principle, we introduce the functional $\Lambda_{t}$ to represent the path-dependent variable linked to the evolution of the spot price $S_{t}$. We then consider evaluation of a derivatives security with the terminal pay-off function $f(S,\Lambda)$. Finally, we extend the joint dynamics in Eq. (6) with the dynamics of augmented variable $\Lambda_{t}=\min_{0\leq t^{\prime}\leq t}S_{t}$:

$$\begin{array}[]{c}S_{t}=rS_{t}dt+\sqrt{V_{t}}S_{t}\left(\rho dB_{t}+\sqrt{1-\rho^{2}}dW_{t}\right),\ \ \ S_{0}=s,\\ dV_{t}=\Phi(V_{t})dt+\Psi(V_{t})dB_{t},\ \ \ V_{0}=v,\\ d\Lambda_{t}=\theta\left(\Lambda_{t}-S_{t}\right)\left(dS_{t}\right)_{-},\ \ \ \Lambda_{0}=S_{0}.\end{array}$$

(35)

The payoff function $f(S,\Lambda)$ is given by:

$$f(S_{T},\Lambda_{T})=\mathbf{1}_{\{\Lambda_{T}>B\}}\Pi\left(S_{T}\right)$$

(36)

The joint dynamic, given by Eqs (35), is Markovian. Accordingly, we can reduce the general formula (34) to the form:

$$P(0,S_{0})=e^{-rT}\int_{V=0}^{\infty}\int_{S=0}^{\infty}\int_{\Lambda=0}^{\infty}f(S,\Lambda)G(T,V,S,\Lambda;0,V_{0},S_{0},\Lambda_{0})d\Lambda dSdV.$$

(37)

Here $G(T,V,S,\Lambda;0,V_{0},S_{0},\Lambda_{0})$ is the risk-neutral probability density function for the joint evolution of state variables in SDE (35). We emphasize that while the payoff function does not depend on the variance $V$, the valuation problem has three spatial variables, including variance, in addition to the time variable.

Below, we introduce a novel method to solve the valuation problem (37) semi-analytically. Using Eq. (22), the density of the price $S_{t}$ is log-normal with time-dependent drift and variance conditional on the filtration generated by stochastic variance $V_{t}$, $\mathbb{F}^{V}$. Therefore, we represent the valuation formula (37) as follows:

$$P(0,S_{0})=e^{-rT}\mathbb{E}_{V}\left[\int_{S=0}^{\infty}\int_{\Lambda=0}^{\infty}f(S,\Lambda)\Gamma\left(\left.T,S,\Lambda;0,S_{0},\Lambda_{0}\right|V(\omega)\right)d\Lambda dS\right],$$

(38)

where $\Gamma$ is the Gaussian density of $S_{t}$ and $\Lambda_{t}$ conditional of the variance path $V(\omega)$, and the expectation is computed over all paths of variance process $V_{t}$. The term in the square brackets is the value of a path-dependent option on $\Lambda$ and $S$, computed in closed-form or numerically.

We apply the MHP to compute the inner integral for barrier options in the closed-form. Thus, we have generalized Eq (27), valid for path-independent options, to path-dependent options, and apply the new result to value barrier options in the semi-closed form.

Eq. (15), conditional on the variance path, can be written as an SDE of the form

$$dX_{t}=\mu\left(t\right)dt+\nu\left(t\right)dW_{t},\ \ \ X_{t}\geq\xi,$$

(39)

where $X_{t}$ the drifted Brownian motion with $\mu,\nu$ given by Eqs (16).

While studying $X_{t}$ on the entire axis $\left(-\infty,\infty\right)$ is almost trivial, dealing with the same process in a semi-bounded domain $\left(\xi,\infty\right)$ might be quite hard. This paper shows how to do it by using the FDM and the MHP.

We can always rescale time $t$ by introducing $\upsilon$, such that $d\upsilon=\nu\left(t\right)dt$. As a result, the governing SDE becomes

$$dX_{\upsilon}=\lambda\left(\upsilon\right)d\upsilon+dW_{\upsilon},$$

(40)

where $\lambda\left(\upsilon\right)=\mu\left(t\right)/\nu\left(t\right)$, $\upsilon\left(t\right)=\int_{0}^{t}v\left(t^{\prime}\right)dt^{\prime}$.

Using the decomposition for $X_{t}$ given by Eq. (20) we present the valuation problems follows:

$$dX_{t}=dY_{t}+dM_{t},\ \ \ Y_{0}=0,\ \ \ M_{0}=0,\ \ \ X_{t}\geq\xi,$$

(41)

where $Y_{t}$ is stochastic part and $M_{t}$ is deterministic part, and $\xi<0$. Then we can express the problem (20) in terms of stochastic variable $Y_{t}$ as follows:

$$dY_{t}=\sqrt{1-\rho^{2}}\sqrt{V_{t}}dW_{t},\ \ \ Y_{0}=0,\ \ \ Y_{t}\geq\xi-M_{t}.$$

(42)

We use a new time variable

$$\upsilon\left(t\right)=\left(1-\rho^{2}\right)I_{t},$$

(43)

and write

$$dY_{\upsilon}=dW_{\upsilon},\ \ \ Y_{0}=0,\ \ \ Y_{\upsilon}\geq\xi-N_{\upsilon},$$

(44)

where

$$N_{\upsilon}=M_{t\left(\upsilon\right)}.$$

(45)

Thus, we have one of the two venues to explore: (A) studying the processes $X_{t}$ or $X_{\upsilon}$ given by Eq. (39) and Eq. (40), respectively; (B) dealing with the processes $Y_{t}$ or $Y_{\upsilon}$ given by Eq. (42) and Eq. (44). In case (A) there is non-zero time-dependent drift and flat boundary; in case (B) there is no drift but the boundary is time-dependent.

Given a variance trajectory, we need to discuss how to calculate $\mu,\nu,\lambda,M$, and $N$. Let $\left\{v_{k}|k=0,1,...,K\right\}$, $v_{0}=v$, be a particular path generated via discretization of Eqs (14) with homogeneous time-step $\Delta t=T/K$. This equation can be discretized in various ways, for example, via the Euler-Maryama scheme or the Milstein scheme. For special cases such as the Feller process corresponding to the Heston model, there are clever schemes tailored to the specific process at hand. However, we are not pursuing them here since it is unnecessary to achieve our objective. We have

$$\begin{array}[]{c}v_{0}=v,\ \ \ v_{k}=v_{k-1}+\kappa\left(\theta-v_{k-1}\right)\Delta t+\varepsilon\sqrt{\Delta t}\eta,\\ I_{0}=0,\ \ \ I_{k}=I_{k-1}+\frac{\Delta t}{2}\left(v_{k}+v_{k-1}\right),\end{array}$$

(46)

$$\begin{array}[]{c}\mu_{k}=r-\frac{\rho\kappa\theta}{\varepsilon}-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)v_{k}+\frac{\rho}{\varepsilon}\frac{\left(\left(v_{k}-v_{k-1}\right)\right)}{\Delta t},\\ \nu_{k}=\sqrt{1-\rho^{2}}\sqrt{v_{k}},\\ M_{0}=0,\ \ \ M_{k}=\frac{kT}{K}\left(r-\frac{\rho\kappa\theta}{\varepsilon}\right)-\left(\frac{1}{2}-\frac{\rho\kappa}{\varepsilon}\right)I_{k}+\frac{\rho}{\varepsilon}\left(v_{k}-v\right).\end{array}$$

(47)

It is clear that

$$\upsilon_{k}=\left(1-\rho^{2}\right)I_{k},\ \ \ \Upsilon=\left(1-\rho^{2}\right)I_{K},$$

(48)

where $\left\{\upsilon_{k}|k=0,1,...,K\right\}$ is an inhomogeneous partition of the interval $\left[0,\Upsilon\right]$. For our purposes, we treat the sequences $\left\{\lambda_{k}=\mu_{k}/\nu_{k}|k=0,1,...,K\right\}$ and
$\left\{M_{k}|k=0,1,...,K\right\}$ as functions of $\upsilon_{k}$, which is possible because $\upsilon_{k}$ is a monotonically increasing sequence, and interpret them accordingly.

We illustrate the corresponding functions in Figures 2, 3.

We emphasize that $\mu$ and $\lambda$ are very irregular, since they depend on the white noise process $dW_{t}/dt$, so that we have to deal with random terms of order $\Delta t^{-1/2}$. At the same time, the moving boundary $M_{t}$ is much more regular, and depends on random terms of order $\Delta t^{1/2}$.

Let us describe how to price a conditional barrier option via the FDM; see Loeper and Pironneau (2009).⁴⁴4Surprisingly, Loeper and Pironneau (2009) do not comment on the irregular nature of the drift coefficient. Specifically, we want to solve the following problem on the semi-axis $\left(\xi,\infty\right)$:

$$\begin{array}[]{c}P_{t}\left(t,X\right)+\mu\left(t\right)P_{X}\left(t,X\right)+\frac{1}{2}\nu\left(t\right)P_{XX}\left(t,X\right)-\varkappa\left(t\right)P\left(t,X\right)=0,\\ P\left(T,X\right)=\Pi\left(X\right),\ \ \ P(t,\xi)=0.\end{array}$$

(49)

We introduce $\bar{P}$, such that

$$P\left(t,X\right)=e^{-\int_{t}^{T}\varkappa\left(t^{\prime}\right)dt^{\prime}}\bar{P}\left(t,X\right),$$

(50)

and get the following problem for $\bar{P}\left(t,X\right)$:

$$\begin{array}[]{c}\bar{P}_{t}\left(t,X\right)+\mu\left(t\right)\bar{P}_{X}\left(t,X\right)+\frac{1}{2}\nu\left(t\right)\bar{P}_{XX}\left(t,X\right)=0,\\ \bar{P}\left(T,X\right)=\Pi\left(X\right),\ \ \ \bar{P}(t,\xi)=0.\end{array}$$

(51)

In the context we are interested in, $\varkappa\left(t\right)=r$, so that $P=\exp\left(r\left(t-T\right)\right)\bar{P}$.

We choose a uniform spatial grid $\left\{\left.x_{m}=\xi+m\left(U-\xi\right)/M\right|m=0,1,...,M\right\}$, where $U$ is a sufficiently remote upper boundary such that $0$ belongs to the spatial grid, $x_{o}=0$, and a temporal grid $\left\{\left.t^{n}\right|t^{0}<t^{1}<...<t^{N}=T\right\}$, which is not necessarily uniform. Then, we apply the usual Crank-Nicolson method for the advection-diffusion diffusion and get the following system of matrix equations for $\bar{P}_{m}^{n}$:

$$\begin{array}[]{c}\bar{P}_{m}^{n+1}-\bar{P}_{m}^{n}+\alpha^{n+1/2}\left(\bar{P}_{m+1}^{n+1}-\bar{P}_{m-1}^{n+1}+\bar{P}_{m+1}^{n}-\bar{P}_{m-1}^{n}\right)\\ +\beta^{n+1/2}\left(\bar{P}_{m+1}^{n+1}-2\bar{P}_{m}^{n+1}+\bar{P}_{m-1}^{n+1}+\bar{P}_{m+1}^{n}-2\bar{P}_{m}^{n}+\bar{P}_{m-1}^{n}\right)=0.\end{array}$$

(52)

where

$$\begin{array}[]{c}\alpha^{n+1/2}=\frac{\mu^{n+1/2}\Delta t^{n+1/2}}{4\Delta x},\ \ \ \beta^{n+1/2}=\frac{\nu^{n+1/2}\Delta t^{n+1/2}}{4\Delta x^{2}}.\end{array}$$

(53)

We emphasize that given the extreme irregularity of $\mu$, using more complicated approaches for treading the drift term is not warranted.

We can write the system of equations in matrix form.

$$\begin{array}[]{c}\mathbb{A}_{\left(\alpha,\beta\right)}^{n+1/2}\bar{P}^{n}=\mathbb{A}_{\left(-\alpha,-\beta\right)}^{n+1/2}\bar{P}^{n+1},\\ \bar{P}_{m}^{N}=\Pi\left(x_{m}\right),\end{array}$$

(54)

where

$$\begin{array}[]{c}\mathbb{A}_{\left(\alpha,\beta,\gamma\right)}=\left(\begin{array}[]{cccccc}1&0&0&0&0&0\\ \alpha-\beta&1+2\beta&-\alpha-\beta&0&0&0\\ .&.&.&.&.&.\\ .&.&.&.&.&.\\ 0&0&0&\alpha-\beta&1+2\beta&-\alpha-\beta\\ 0&0&\beta&-\alpha-4\beta&4\alpha+5\beta&1-3\alpha-2\beta\end{array}\right).\end{array}$$

(55)

For brevity, the superscripts are omitted. At infinity we apply one-sided derivatives and use the PDE itself to formulate the boundary condition.

Alternatively, we can rescale time, reduce the problem (29) to the form:

$$\begin{array}[]{c}\bar{P}_{\upsilon}+\lambda\left(\upsilon\right)\bar{P}_{x}+\frac{1}{2}\bar{P}_{xx}=0,\\ \bar{P}\left(\Upsilon,x\right)=\Pi\left(x\right),\ \ \ \bar{P}(\upsilon,\xi)=0,\end{array}$$

(56)

and apply the Crank-Nicolson method to problem (56).

Of course, we can attack the pricing problem by solving the corresponding Fokker-Planck equation for Green’s function $G$:

$$\begin{array}[]{c}G_{t}(t,X)+\mu\left(t\right)G_{X}(t,X)-\frac{1}{2}\nu\left(t\right)G_{XX}(t,X)+\varkappa\left(t\right)G(t,X)=0,\\ G\left(0,X\right)=\delta\left(X\right),\ \ \ G(t,\xi)=0,\end{array}$$

(57)

and write

$$\begin{array}[]{c}P\left(0,0\right)=\int\limits_{0}^{\infty}G\left(0,0;T,x\right)\Pi\left(x\right)dx.\end{array}$$

(58)

We introduce $\bar{G}$, such that

$$G\left(t,X\right)=e^{-\int_{0}^{t}\varkappa\left(t^{\prime}\right)dt^{\prime}}\bar{G}\left(t,X\right).$$

(59)

Accordingly,

$$\begin{array}[]{c}\bar{G}_{t}(t,X)+\mu\left(t\right)\bar{G}_{X}(t,X)-\frac{1}{2}\nu\left(t\right)\bar{G}_{XX}(t,X)=0,\\ \bar{G}\left(0,X\right)=\delta\left(X\right),\ \ \ \bar{G}(t,\xi)=0,\end{array}$$

(60)

$$\begin{array}[]{c}P\left(0,0\right)=e^{-\int_{0}^{T}\varkappa\left(t^{\prime}\right)dt^{\prime}}\int\limits_{0}^{\infty}\bar{G}\left(0,0;T,x\right)\Pi\left(x\right)dx.\end{array}$$

(61)

As before, we apply the Crank-Nicolson method for the advection-diffusion and get the following system of matrix equations for $\bar{G}_{m}^{n}$:

$$\begin{array}[]{c}\bar{G}_{m}^{n+1}-\bar{G}_{m}^{n}+\alpha^{n+1/2}\left(\bar{G}_{m+1}^{n+1}-\bar{G}_{m-1}^{n+1}+\bar{G}_{m+1}^{n}-\bar{G}_{m-1}^{n}\right)\\ -\beta^{n+1/2}\left(\bar{G}_{m+1}^{n+1}-2\bar{G}_{m}^{n+1}+\bar{G}_{m-1}^{n+1}+\bar{G}_{m+1}^{n}-2\bar{G}_{m}^{n}+\bar{G}_{m-1}^{n}\right)=0.\end{array}$$

(62)

In the matrix form we get

$$\begin{array}[]{c}\mathbb{A}_{\left(\alpha,\beta\right)}^{n+1/2}\bar{G}^{n+1}=\mathbb{A}_{\left(-\alpha,-\beta\right)}^{n+1/2}\bar{G}^{n},\\ \bar{G}^{0}=\delta_{o,m}.\end{array}$$

(63)

Once $\bar{G}^{N}$ is found, we can represent $P_{k}^{0}$ as

$$\begin{array}[]{c}P_{k}^{0}=e^{-rT}\Delta x\sum\limits_{m=1}^{M}\bar{G}_{m}^{N}\Pi_{m}.\end{array}$$

(64)