In the previous sections we have only considered price processes where the coefficients in the governing SDE are constant (see section on jump-diffusion models). One can also consider models where the coefficients are functions of time and the stock price.
An alternative approach is to allow the paremeters such as the volatility \(\sigma\) to be random variables themselves. More precisely, we look at models of the form \((Y_t, V_t)\) where the log-price \(Y_t\) and the volatility \(V_t\) are governed by different SDEs. These type of models are often called stochastic volatility models and an example of which is the Barndorff-Nielson--Shephard (BNS) model. The BNS model is of the form

\(\begin{align} \mathrm{d}Y_t &= \Big(r -\delta - \frac{1}{2}V_t\Big)\mathrm{d}t + \sqrt{V_t}\mathrm{d}W_t + \rho\mathrm{d}J_{\lambda t}, \\ \mathrm{d}V_t &= -\lambda V_t \mathrm{d}t + \mathrm{d}J_{\lambda t}, \end{align}\)

where \(r,\lambda>0\), \(\rho<0\), \((W_t)_{t\geq 0}\) is a standard Brownian motion and \((J_{\lambda t})_{t\geq 0}\) is a pure jump process that increases almost surely and has Levy measure \(\nu\). \(\delta\) is a time-independent constant which is chosen so that the discounted asset price \(\tilde{S}_t = \exp(-rt + Y_t)\) is a martingale. It can be shown that necessarily \(\delta = \lambda\kappa(\rho)\) where \(\kappa\) is the cumulant generating function of \((J_t)_{t\geq 0}\) given by

\( \kappa(u) = \int_0^\infty (e^{u\xi} - 1) \;\nu(\mathrm{d}\xi). \)

The negative constant \(\rho\) is there to introduce a correlation between the log-price and the volatility; a jump in the volatility will cause a jump downwards in the log-price. Typically, the jump term \((J_t)_{t\geq 0}\) is taken to be a compound Poisson process

\( J_t = \sum_{i=1}^{N_t} U_i \)

where \( (N_t)_{t\geq 0} \) is a Poisson process with intensity \(\alpha\) and \(\{U_i\}_{i\in\mathbb{N}}\) is a sequence of i.i.d. random variables which we assume to be exponentially distributed with mean \(1/\eta\) and so \(J_{\lambda t}\) has Levy measure given by \(\nu(\mathrm{d}h) = \alpha\lambda g(h)\) with \(g(h) = \eta\exp(-\eta h)\) being the density of a exponentially distributed random variable. For this choice of \( (J_t)_{t\geq 0}\) a simple integration shows that

\(\kappa(\rho) = \frac{\alpha \rho}{\eta - \rho}\).

As before option prices are given as expectations of the payoff under the risk neutral measure so Monte Carlo methods can be used.Notice that the BNS model is a two-dimensional model so one needs to approximate it by a Markov chain on a finite two-dimensional state space. In one-dimensional models the Q-matrix \(\Lambda\) of the approximating chain is a \(N\times N\) matrix where \(N\) is the number of grid points. If now we have \(N\) grid points in bothe the log-price and the volatility component then the Q-matrix is a \(N^2\times N^2\) matrix. Thus, simulations of these higher dimensional models can be computationally expensive.

As in the jump-diffusion case we calculate the Q-matrix of the Markov chain by considering the jump and the diffusion components of the process separate and compute the corresponding transition rate matrices \(\Lambda_J\) and \(\Lambda_D\), then \(\Lambda\) is defined as the sum \(\Lambda\) = \(\Lambda_D + \Lambda_J\).

Diffusion Matrix

We can express the diffusion part of the BNS model in vector form as

\(\mathrm{d}\begin{pmatrix} Y_t \\ V_t \end{pmatrix} = \begin{pmatrix} r- \delta - \frac{1}{2}V_t \\ -\lambda V_t \end{pmatrix}\mathrm{d}t + \begin{pmatrix} \sqrt{V_t} \\ 0 \end{pmatrix}\mathrm{d}W_t\)

which provides the infinitesimal generator \(\mathcal{L}_D\) as

\(\mathcal{L}_Df(y,v) = \Big(r - \delta - \frac{1}{2}v\Big)\frac{\partial f}{\partial y}(y,v) - \lambda v\frac{\partial f}{\partial v}(y,v) + \frac{1}{2}v\frac{\partial^2 f}{\partial y^2}(y,v).\)

where \(f \in C_0^2 (\mathbb R ^2)\) the space of twice differentiable functions such that \(\lim_{|x|\to\infty}f(x) = 0 \). One can discretize this PDE by finite differences using either forward or backward difference for the first order derivatives:

\(\frac{1}{h}\Big( f(y + h,v) - f(y,v)\Big)\quad\text{or}\quad \frac{1}{h}\Big( f(y,v) - f(y - h, v)\Big)\)

The choice of which one to use is determined by the condition that the result must be a Q-matrix. For the second order derivative term central difference can be used:

\(\frac{1}{h^2}\Big( f(y + h, v) - 2f(y,v) + f(y - h,v)\Big)\).

The forward, backward and central differences in the \(v\) variable are defined similarly. Let \(\mathbb{G} = \{(y_i,v_j) : 1\leq i,j\leq N, y_1 > \cdots > y_N, v_1 < \cdots < v_N\}\) be a grid with grid spacing \(h_v\) in the volatility component and grid spacing \(h_y = |\rho|h_v\) in the log-price component. The choice of the grid spacing \(h_y\) is so that a jump in the volatility of one grid point translates to a jump of one grid point in the log-price and so on. Now denote \(f(y_i,v_j)\) by \(f_{ij}\) then the discretisation of \(\Lambda_D\) is

\(\frac{1}{h_y}\Big(\mu f_{i-1j} - \mu f_{ij} - \frac{1}{2}v_jf_{ij} + \frac{1}{2}v_jf_{i+1}j\Big) + \frac{1}{h_v}\Big(\lambda v_jf_{ij-1} - \lambda v_jf_{ij}\Big) + \frac{1}{h_y^2}\Big(\frac{1}{2}v_jf_{i-1j} - v_jf_{ij} + \frac{1}{2}v_jf_{i+1j}\Big)\)

It is easy to see that the above equation can be expressed as a matrix-vector product \(\Lambda_J \mathbf{f}\) where \( \mathbf{f} = \big(f_{11},f_{12},\ldots,f_{1N},f_{21},\ldots,f_{NN}\big)^T \) and the \(\Lambda_J\) is the desired discretisation.

Jump Matrix

In general, the generator of a compound Poisson process on \(\mathbb{R}^d\) is given by

\( \mathcal{L}_J f(x) = \int_{\mathbb{R}^d} \big(f(x+h) - f(x)\big)\;\nu(\mathrm{d}h) \),

where \(f \in C_0(\mathbb{R}^d)\). In the BNS model, since the log-price process and the volatility are both driven by the same jump process with correlation \(\rho < 0\) we can write the generator \(\mathcal{L}_J\) for the jump component as

\(\mathcal{L}_Jf(y,v) = \int_{\mathbb{R}} \big( f(y + \rho h, v + h) - f(y,v) \big)\;\nu(\mathrm{d}h),\)

The above integral can be approximated by Riemann sums on the grid \(\mathbb{G}\),

\(\begin{align*}\mathcal{L}_J f(y_i,v_j) &\approx \sum_{n=1}^{N - (i\vee j)} \Big( f(y_i + n\rho h_v,v_j + nh_v) - f(y_i,v_j) \Big) \alpha\lambda g(h_v)nh_v \\ &= \sum_{n=1}^{N - (i\vee j)} \Big( f(y_{i + n} ,v_{j + n}) - f(y_i,v_j) \Big) \alpha\lambda g(h_v)nh_v \end{align*}\).

Again the above can be written as a matrix-vector product \(\Lambda_J\mathbf{f}\) where \(\mathbf{f}\) is as above.

 

The above graphs show a sample path of the log-price and the volatility obtained from the approximating Markov chain. A grid of size 100 \(\times\) 100 was used to generate the path. Notice that the log-price looks like a typical jump-diffusion process while the volatility maoves up entirely by jumps and slowly drifts to zero in the abscence of jumps due to the negative term in the governing SDE.