Multinomial method for option pricing under Variance Gamma

The VG process is obtained by time changing a Brownian motion with drift. The new time variable is a stochastic process $T_{t}\sim\Gamma(\mu t,\kappa t)$ with density²²2Usually the Gamma distribution is parametrized by a shape and scale positive parameters $T\sim\Gamma(\rho,\zeta)$. The Gamma process $T_{t}\sim\Gamma(\rho t,\zeta)$ has pdf $f_{T_{t}}(x)=\frac{\zeta^{-\rho t}}{\Gamma(\rho t)}x^{\rho t-1}e^{-\frac{x}{\zeta}}$ and has moments $\mathbb{E}[T_{t}]=\rho\zeta t$ and $\mbox{Var}[T_{t}]=\rho\zeta^{2}t$. Here we use a parametrization as in [18] such that $\mathbb{E}[T_{t}]=\mu t$ and $\mbox{Var}[T_{t}]=\kappa t$, so $\zeta=\frac{\kappa}{\mu}$, $\rho=\frac{\mu^{2}}{\kappa}$.:

The Gamma process $T_{t}$ is a subordinator. A subordinator is a one dimensional Lévy process that is non-decreasing almost surely. Therefore it is suitable to represent a time variable. It is possible to prove that every subordinator is a finite variation process (see [4]).

Considering a Brownian motion with drift $X_{t}=\theta t+\sigma W_{t}$, with $W_{t}\sim\mathcal{N}(0,t)$, let’s replace the time variable by the Gamma subordinator $T_{t}\sim\Gamma(t,\kappa t)$ (with $\mu=1$). We obtain the Variance Gamma process:

It depends on three parameters:

• •

  $\theta$, the drift of the Brownian motion,

• •

  $\sigma$, the volatility of the Brownian motion,

• •

  $\kappa$, the variance of the Gamma process.

The probability density function of $X_{t}$ can be computed conditioning on the realization of $T_{t}$:

where the function $K$ is a modified Bessel function of the second kind (see [18] for the computations). The characteristic function can be obtained by conditioning too:

see proposition 1.3.27 in [4]. The Lévy symbol is:

Using the formula (70) in Appendix B for the cumulants we derive:

The VG Lévy measure is³³3In [18] the authors derive the expression for the Lévy measure by representing the VG process as the difference of two Gamma processes.:

and it satisfies conditions (4) and (5). The VG process can be represented as a pure jump process as in (6) and (7), with no additional drift $b^{\prime}=0$.

All the informations are contained in the Lévy measure (15), which completely describes the process. Even if the process has been created by Brownian subordination, it has no diffusion component. The Lévy triplet is $(\int_{|x|<1}x\nu(dx),0,\nu)$. Using the formalism of Poisson integrals in Appendix A, the Lévy symbol (13) has the representation⁴⁴4See Example 8.10 in [22].:

Under the risk neutral measure $\mathbb{Q}$, the dynamics of the stock price is described by an exponential Lévy model:

where $r$ is the risk free interest rate, and $L_{t}$ is a general Lévy process. Under $\mathbb{Q}$, the discounted price is a $\mathbb{Q}$-martingale:

and thus $\mathbb{E}^{\mathbb{Q}}[e^{L_{t}}|L_{0}=0]=1$. The condition for the existence of the exponential moment $\mathbb{E}[e^{L_{t}}]<\infty$ is equivalent to:

as proved in Lemma 25.7 in [22]. The VG process $X_{t}$ has finite exponential moment. In order to satisfy the martingale condition⁵⁵5 To obtain the correction term $\omega$ we have to find the exponential moment of $X_{t}$ using its characteristic function: $\displaystyle\mathbb{E}[e^{X_{t}}]=\phi_{X_{t}}(-i)=e^{-\omega t}$ . we need to add a correction term to $X_{t}$. The following process is a martingale:

where $w=\frac{1}{\kappa}\log(1-\theta\kappa-\frac{1}{2}\sigma^{2}\kappa)$. Passing to the log-price $Y_{t}=\log(S_{t})$, we get a process in the form of Eq. (6) with $b^{\prime}=r+\omega$:

Let $V(t,Y_{t})$ be the value of an option at time $t$. We assume that $V(t,y)\in C^{1,1}([t_{0},T],\mathbb{R})$ and has a polynomial growth at infinity.
By the martingale pricing theory, the discounted price of the option is a martingale and it is possible to derive the PIDE for the price of the option:

where $\mathcal{L}^{Y_{t}}$ is the infinitesimal generator of the log-price process (22). The resulting PIDE is:

To determine the parameters of the Markov chain we require that its local moments are equal to that of the continuous process. First, we rewrite the continuous process (22) as the sum of a drift term and a martingale term:

where $\theta=\int_{\mathbb{R}\backslash\{0\}}x\nu(dx)=\mathbb{E}\bigl{[}\int_{\mathbb{R}\backslash\{0\}}xN(1,dx)\bigr{]}$ is the expected value of the VG process in 16, when $\Delta t=1$. The integral with respect to the compensated Poisson measure $\tilde{N}(\Delta t,dx)$ is a martingale (see Appendix A).

We can pass to log-prices $Y_{n}=\log(S_{n})$ in the discrete Eq. (25), and write it as the sum of a drift component and a random variable with $L$ possible outcomes:

The term $\bar{b}\,\Delta t$ is the drift term, while $l$ is a random variable that assumes values in $[1,2,...,L]$ with probability $p_{l}$. It has to satisfy the martingale condition:

with $\alpha(\Delta t)$ a function of $\Delta t$.

The corresponding up/down factors have the following representation:

and we can readily see that if $u/d$ is constant the tree recombines.

Given the mean $c_{1}=\mathbb{E}[\Delta Y]=\bar{b}\Delta t$, the $k$-central moment is:

The moment matching condition requires that the central moments of the discrete process (28) are equal to the central moments of the continuous process (27):

We fix $L=5$, and using the relation between central moments and cumulants (Eq. (71) in Appendix B) we solve the linear system of equations for the transition probabilities:

The drift parameter is $\bar{b}=r+\omega+\theta$. The only missing parameter to determine is $\alpha(\Delta t)$. This is a function of the time increment $\Delta t$ and can be determined using the higher order terms in the moment matching condition together with the condition of positive probabilities.

Recall that the well known binomial model [10] assumes the value $\alpha(\Delta t)=\tilde{\sigma}\sqrt{\Delta t}$, that represents the volatility of the increments in the time interval $\Delta t$. In the trinomial model, the parameter $\alpha(\Delta t)$ assumes value $\alpha(\Delta t)=\frac{1}{2}\tilde{\sigma}\sqrt{3\Delta t}$, see for instance [25]. For the multinomial method a good representation for $\alpha(\Delta t)$ is:

where $\bar{\kappa}=c_{4}/c_{2}^{2}$ is the excess of kurtosis⁶⁶6We use the bar over $\kappa$, to distinguish the kurtosis from the variance of the gamma process $\kappa$.. We refer to [27] for the derivation. This choice guarantees that the probabilities $p_{i}$ for $i=1...5$ are always positive and sum to one. We can replace the expression (33) inside (32), to obtain the simpler form:

where $s=c_{3}/\sqrt{c_{2}^{3}}$ is the skewness.

We call a generic jump process (6) with first four cumulants $c_{1}$,$c_{2}$,$c_{3}$,$c_{4}$ equal to the VG cumulants (14), the approximated process $X^{A}$. The cumulant generating function of the increment $\Delta X^{A}$ has the following series representation (see Appendix (B)):

We are interested in the approximation of a VG process with drift (27), therefore we require that $c_{1}=\bar{b}\Delta t=(r+\omega+\theta)\Delta t$.

For the proof of this theorem we refer to Section 4.2 of [27]. The authors first define the jump sizes and their respective probabilities, and then prove that when $\Delta t\to 0$ the characteristic function of the pentanomial tree converges to the characteristic function of the compound Poisson process.

We suggest the following algorithm for pricing with the multinomial method:

1. 1.

   Compute the transition probabilities vector (34).

2. 2.

   Compute the up/down factors $u$ and $d$ (29) and the vector of prices $S_{N}$ at terminal time $N$ as in Eq. (26).

3. 3.

   Evaluate the payoff of the option $V^{N}(S_{N})$ at terminal time $N$.

4. 4.

   Given the option’s values at time $t_{n+1}$ compute the values at time $t_{n}$. The value is the conditional expectation:

   $$V^{n}(s^{(k)}_{n})=e^{-r\Delta t}\mathbb{E}^{\mathbb{Q}}\biggl{[}V^{n+1}(S_{n+1})\bigg{|}S^{(k)}_{n}=s^{(k)}_{n}\biggr{]}.$$

   (38)

5. 5.

   If computing the price of an American option, the value at the previous time level is the maximum between the conditional expectation and the intrinsic value of the option. For an American put we have:

   $$V^{n}(s^{(k)}_{n})=\max\biggr{\{}e^{-r\Delta t}\mathbb{E}^{\mathbb{Q}}\biggl{[}V^{n+1}(S_{n+1})\bigg{|}S^{(k)}_{n}=s^{(k)}_{n}\biggr{]},K-s^{(k)}_{n}\biggr{\}}.$$

   (39)

6. 6.

   Iterate the algorithm until the initial time $t_{0}$.

In the parameters calibration procedure, sometimes it is common to estimate first the historical parameters, and use them as initial guess for the least squares minimization that recovers the risk neutral parameters. In [23] are presented several methods for historical parameters estimation of the VG density. We use the simple method of moments to estimate the parameters for the data in Fig. 1. In all future calculations we consider the risk neutral VG parameters in Table 4.1.

We compare the numerical results for European call and put options obtained with the multinomial and the PIDE approaches.

• •

  VG PIDE: We solve the VG PIDE following the method proposed by [9]. The Lévy measure is singular in the origin and this is a problem for the computation of the integral term in (24). Fixing a truncation parameter $\epsilon>0$, we approximate the infinite activity martingale jump component with sizes smaller than $\epsilon$, with a Brownian motion with same variance. The resulting approximated PIDE has the “jump-diffusion-like” form:

  		$\displaystyle\frac{\partial V(t,x)}{\partial t}+\bigl{(}r-\frac{1}{2}\sigma_{\epsilon}^{2}-w_{\epsilon}\bigr{)}\frac{\partial V(t,x)}{\partial x}+\frac{1}{2}\sigma_{\epsilon}^{2}\frac{\partial^{2}V(t,x)}{\partial x^{2}}$		(40)
  		$\displaystyle+\int_{|z|\geq\epsilon}V(t,x+z)\nu(dz)=(\lambda_{\epsilon}+r)V(t,x).$

  where we introduced the parameters $\sigma_{\epsilon}^{2}=\int_{|z|<\epsilon}z^{2}\nu(dz)$, $\omega_{\epsilon}=\int_{|z|\geq\epsilon}(e^{z}-1)\nu(dz)$ and $\lambda_{\epsilon}=\int_{|z|\geq\epsilon}\nu(dz)$. More details are in [8]. We solve the PIDE (40) using the implicit-explicit finite difference scheme proposed in [9], and choosing the truncator parameter $\epsilon=1.5\delta x$, where $\delta x$ is the size of the space step. It turns out that the solution of the discretized equation convergences very slowly to the option price, and therefore we required a grid with $14000$ space steps and $7000$ time steps. The algorithm is written in Matlab and runs on an Intel i7 (7th Gen) with Linux. It takes about 30 minutes to complete.

• •

  Multinomial: We follow the algorithm proposed in the previous section. The number of time steps for all the computations is $N=2000$. In the table 4.2 we show a convergence table for the prices of European calls, puts and American puts. It is straightforward to see that the convergence is quite fast.

Figures (2) and (3) show the prices obtained by the multinomial method compared with those obtained by PIDE. In table 4.2 we compare directly the call/put numerical values obtained with the two methods.

Pricing vanilla call and put European options is quite simple and the best approach is to use the closed formula derived in [18]. The big advantage of the multinomial method is in the computation of American options prices, where there is no closed formula and all the other approaches, such as PIDEs and Least Squares Monte Carlo, are difficult to implement and much slower.

In this section we present the numerical results obtained with the multinomial method algorithm for American put options, and compare them with the PIDE method (see fig. 4). The PIDE (40) is modified in order to take in account the early exercise feature:

To solve this equation we use the same settings used for the Eq. (40): an implicit-explicit finite difference scheme, with 14000 space steps and 7000 time steps. The algorithm takes about 33 minutes to run.

The numerical values of the prices obtained with the multinomial and PIDE methods are collected in Tab. 4.3. The run times for the multinomial algorithm are shown in the convergence table 4.2.

In Table 4.3, we consider also European and American put option prices calculated with the Black-Scholes (BS) models. The BS volatility is chosen equal to the VG volatility $\sigma^{BS}=(\sigma^{2}+\theta^{2}\kappa)$. As expected, the deep OTM (out of the money) prices computed under VG are higher than the corresponding prices computed under BS. This is a consequence of the shape of the VG density function (10), which has heavier tails than the normal distribution. This means that the probability of a deep OTM option to return in the money, is higher if calculated with the VG model than BS, and therefore we get higher option prices.

The Black-Scholes prices are computed using a binomial algorithm. For European options the prices converge to the BS closed formula prices. The same values can be obtained using the multinomial algorithm for the VG process and setting $\theta=\kappa=0$ and $\sigma=\sigma^{BS}$. Recall that under the Black-Scholes model, the log-returns follow a Brownian motion. Looking at the definition of the VG process (9), it is evident that when $\theta$ and $\kappa$ are zero, the process becomes a Brownian motion:

As a consequence, the price process (21) converges to the Geometric Brownian Motion:

where: