Variance Gamma Model: Characteristic Function and Python

In October 1990 Dilip Madan and Eugene Seneta published a return model with no Brownian component in calendar time. Eight years later Madan, Peter Carr and Eric Chang turned it into the option-pricing model in The Variance Gamma Process and Option Pricing.

Variance gamma, or VG, is built by running a Brownian motion on a random gamma clock. The construction produces skew and excess kurtosis with only three parameters, while retaining a closed-form density and characteristic function.

It sits between the finite-activity Merton jump diffusion and more flexible Lévy families. The Carr-Madan FFT turns its characteristic function into a fast vanilla pricer.

Gamma time change

Let $G_t$ be a gamma process normalised by

$$G_t\sim\operatorname{Gamma}\!\left(\frac{t}{\nu},\nu\right), \qquad \mathbb E[G_t]=t, \qquad \operatorname{Var}(G_t)=\nu t,$$

where the second gamma argument is the scale. Let $W$ be an independent Brownian motion and define

$$X_t=\theta G_t+\sigma W_{G_t}.$$

Conditional on $G_t=g$, $X_t$ is normal with mean $\theta g$ and variance $\sigma^2g$. Integrating over the gamma clock gives the VG process.

The clock supplies useful intuition. When $G_t$ exceeds calendar time, the Brownian driver experiences more business time and produces a more dispersed return. When it is small, the return is compressed. The parameter $\nu$ controls variability of this clock because $\operatorname{Var}(G_t)=\nu t$. As $\nu\to0$, $G_t$ concentrates at $t$ and

$$X_t\longrightarrow\theta t+\sigma W_t.$$

For positive $\nu$, mixing normals with different conditional variances creates heavy tails. The drift $\theta G_t$ correlates clock realisations with return direction: negative $\theta$ makes high-activity periods more negative on average and creates left skew. This is more precise than saying that $\theta$ is simply a skew parameter, since it also changes the mean and variance.

It is a Lévy process with stationary independent increments. In calendar time it is pure jump, has infinitely many small jumps on every interval and has finite variation. The Brownian motion in the construction does not leave a continuous Brownian component after gamma subordination.

If time is measured in years, $\nu$ has units of years, $\theta$ has units per year and $\sigma$ has units per square root year. $X_t$ is a dimensionless log return.

Characteristic function

Conditioning on $G_t$ gives

$$\begin{aligned} \phi_X(u;t) &=\mathbb E[e^{iuX_t}]\\ &=\mathbb E\!\left[ e^{(iu\theta-\frac12\sigma^2u^2)G_t} \right]\\ &=\left( 1-i\theta\nu u+\frac12\sigma^2\nu u^2 \right)^{-t/\nu}. \end{aligned}$$

This convention places $\theta$ before $\nu$ in the imaginary term. Changing the Fourier sign convention changes that sign too; copying a characteristic function without its convention is a common source of reflected skew.

The conditioning step also supplies an exact simulator for European payoffs. Draw $G_T$ from a gamma distribution with shape $T/\nu$ and scale $\nu$, draw an independent standard normal $Z$, then set

$$X_T=\theta G_T+\sigma\sqrt{G_T}Z.$$

No calendar-time discretisation is needed. This makes Monte Carlo a useful independent check on Fourier prices, although it converges too slowly to certify tiny out-of-the-money values without variance reduction.

The cumulants are more informative than loose parameter labels:

$$\mathbb E[X_t]=\theta t, \qquad \operatorname{Var}(X_t)=(\sigma^2+\theta^2\nu)t,$$ $$\operatorname{Skew}(X_t)= \frac{(3\sigma^2\theta\nu+2\theta^3\nu^2)t} {[(\sigma^2+\theta^2\nu)t]^{3/2}},$$

and

$$\operatorname{ExKurt}(X_t)= \frac{3\sigma^4\nu+12\sigma^2\theta^2\nu^2+6\theta^4\nu^3} {(\sigma^2+\theta^2\nu)^2t}.$$

Thus $\sigma$ is not the annual standard deviation unless $\theta=0$ or $\nu\to0$. Nor do the parameters act independently: changing $\theta$ also changes variance and kurtosis.

Martingale correction

Suppose the three VG parameters are already risk-neutral. Set

$$S_t=S_0\exp\!\left((r-q+\omega)t+X_t\right).$$

Since

$$\mathbb E[e^{X_t}] =\left(1-\theta\nu-\frac12\sigma^2\nu\right)^{-t/\nu},$$

the discounted total-return stock is a martingale when

$$\boxed{ \omega=\frac1\nu \log\!\left(1-\theta\nu-\frac12\sigma^2\nu\right) }.$$

The exponential moment exists only if

$$1-\theta\nu-\frac12\sigma^2\nu>0.$$

This is a real restriction, not something that “always holds” in calibration. The log-price characteristic function is

$$\phi_{\log S_t}(u)= e^{iu[\log S_0+(r-q+\omega)t]} \left(1-i\theta\nu u+\frac12\sigma^2\nu u^2\right)^{-t/\nu}.$$

More generally,

$$\mathbb E[e^{pX_t}] =\left(1-\theta\nu p-\frac12\sigma^2\nu p^2\right)^{-t/\nu},$$

where the quadratic inside brackets is positive. This moment interval matters twice: $p=1$ is needed for the stock martingale, while Carr-Madan damping needs $p=\alpha+1$. For the numerical parameters used below, the quadratic roots are $-12.5$ and $20$. Thus $\alpha=1.5$ is safely admissible and positive damping requires $0<\alpha<19$. This check should be performed from calibrated parameters for every maturity, not assumed from a default damping value.

Like other exponential Lévy models, VG does not have a unique equivalent martingale measure. An Esscher transform, a minimal-entropy choice or direct risk-neutral calibration can produce different pricing parameters. The martingale correction fixes drift after that choice; it does not choose the measure.

Verified FFT code

The Simpson weights below use zero-based indexing. Their first coefficients are $1/3,4/3,2/3$, not $1,2/3,4/3$.

import numpy as np

def vg_cf(u, S0, T, r, q, sigma, nu, theta):
    base = 1 - theta*nu - 0.5*sigma**2*nu
    if base <= 0:
        raise ValueError("The martingale exponential moment does not exist")
    omega = np.log(base) / nu
    vg = (1 - 1j*theta*nu*u + 0.5*sigma**2*nu*u**2)**(-T/nu)
    return np.exp(1j*u*(np.log(S0) + (r - q + omega)*T)) * vg

def vg_call_grid(S0, T, r, q, sigma, nu, theta,
                 alpha=1.5, N=4096, eta=0.25):
    p = alpha + 1
    if 1 - theta*nu*p - 0.5*sigma**2*nu*p**2 <= 0:
        raise ValueError("The damping moment E[S**(alpha+1)] does not exist")
    dk, j = 2*np.pi/(N*eta), np.arange(N)
    b, v = 0.5*N*dk, eta*j
    k = -b + dk*j
    u = v - 1j*(alpha + 1)
    psi = (np.exp(-r*T) * vg_cf(u, S0, T, r, q, sigma, nu, theta)
            / (alpha**2 + alpha - v**2 + 1j*(2*alpha + 1)*v))
    weight = (3.0 - (-1.0)**j - (j == 0)) / 3.0
    values = np.fft.fft(np.exp(1j*b*v) * psi * eta * weight)
    calls = np.exp(-alpha*k) * values.real / np.pi
    return k, calls

The function returns log-strikes. Use, for example, np.interp(np.log(K), k, calls) for an on-grid-range strike, then repeat with other $N$, $\eta$ and admissible $\alpha$ values.

With $S_0=K=100$, $T=0.5$, $r=3\%$, $q=0$, $\sigma=0.20$, $\nu=0.20$ and $\theta=-0.15$, we obtain

$$\omega=0.1283387, \qquad \sqrt{\sigma^2+\theta^2\nu}=21.095\%.$$

The interpolated VG call is $6.36371$. Black-Scholes at $20\%$ is $6.37103$, while Black-Scholes at the VG process standard deviation is $6.67519$. Fat tails do not mechanically raise every option price; the martingale correction and negative skew redistribute probability.

Calibration and checks

For one maturity, fit $(\sigma,\nu,\theta)$ to liquid option prices or implied volatilities. Price errors, volatility errors and relative errors are not interchangeable. If fitting prices, normalising by Black-Scholes vega approximates a volatility-space objective; inverse bid-ask variance reflects quote uncertainty.

A single constant parameter triple determines every maturity. Refitting $\nu$ separately at each expiry may improve slices but is no longer one VG Lévy process, because stationary increments have been abandoned. Decide whether the goal is a slice interpolator or a dynamic model before calibrating.

Useful tests include:

• $\phi(0)=1$ and $\phi(-i)=S_0e^{(r-q)T}$;
• the limit $\nu\to0$, which gives Brownian log returns with drift $\theta$ before martingale correction;
• call bounds, strike monotonicity, convexity and put-call parity;
• agreement between FFT, direct quadrature and the closed-form VG formula in Madan, Carr and Chang;
• stability across transform grids and damping values inside the finite moment strip.

The Black-Scholes model is a useful limiting and numerical benchmark.

Calibration failure modes and practitioner choices

The three parameters are parsimonious but not orthogonal. Near $\nu=0$, many combinations produce distributions close to Brownian motion, so $\nu$ can be weakly identified from a narrow strike range. A more negative $\theta$ changes skew, variance and the martingale correction simultaneously. Increasing $\sigma$ can offset part of the central curvature created by $\nu$. Optimisers may therefore find several shallow minima with similar prices.

Use parameter transforms or explicit bounds to enforce $\sigma>0$, $\nu>0$ and the exponential-moment condition. A generic positivity constraint on each parameter is wrong because $\theta$ is signed. Reject trial points outside the required moment strip before evaluating complex powers. Letting invalid points return large finite prices can create discontinuous objectives and misleading finite-difference gradients.

Calibration weights should reflect the task. Price-space least squares favours expensive options. Implied-volatility errors become unstable where vega is tiny. Inverse bid-ask variance is defensible when spreads represent observation uncertainty, while vega-normalised price errors approximately target volatility fit in liquid regions. Whatever objective is chosen, report residuals in price, implied volatility and fractions of bid-ask spread.

Term structure is the central modelling decision. A single Lévy triplet across maturities preserves stationary independent increments but often misses the changing equity skew. Independent parameter triples by expiry can fit better, yet they define a collection of marginal distributions rather than one VG process. Before using such slices for path-dependent products, specify an interpolation or time-inhomogeneous process and retest calendar arbitrage. A good static fit does not supply coherent dynamics by itself.

A thorough numerical audit should include:

• gamma draws with empirical mean $T$ and variance $\nu T$;
• simulated $X_T$ moments against the analytic mean and variance;
• the martingale identity $\phi_{\log S_T}(-i)=S_0e^{(r-q)T}$;
• agreement among FFT, direct quadrature and exact-terminal Monte Carlo;
• call bounds, strike monotonicity, convexity and put-call parity;
• stability as $N$, $\eta$ and admissible $\alpha$ change;
• convergence to Black-Scholes as $\nu\to0$, with volatility and drift interpreted after applying $\omega$.

The limiting test deserves care. As $\nu\to0$,

$$\omega\to-\theta-\frac12\sigma^2,$$

so $(r-q+\omega)t+X_t$ tends to $(r-q-\sigma^2/2)t+\sigma W_t$. The $\theta t$ terms cancel. A test that expects a residual stock drift of $\theta$ after martingale correction is testing the wrong limit.

For risk management, inspect parameter profiles and day-to-day stability, not only the final objective value. If $\nu$ jumps while fitted vanillas barely move, tail-sensitive exotic values can still change materially. Stress the moment-strip boundary too. Calibration close to it can make high-order moments and damped transforms unstable even while ordinary call prices appear reasonable.

What VG misses

Independent increments rule out volatility clustering. A large move today does not make a large move tomorrow more likely. Constant parameters also struggle with the joint term structure of index smiles.

VG's infinite activity describes an accumulation of small jumps, while Merton's finite-activity process is a cleaner model of isolated news events. Neither distinction guarantees a better market fit. CGMY adds a fourth parameter for jump activity and separate tail control, but brings further identification risk.

Stochastic clocks can add persistence to a Lévy driver. That improves dynamics at the cost of the three-parameter simplicity that makes VG attractive in the first place.

Frequently Asked Questions

Is variance gamma a pure-jump process?

Yes. It has no continuous Brownian component in calendar time, despite being constructed from a Brownian motion. It has infinite jump activity and finite variation.

What is the variance of a VG return?

For horizon $t$, it is $(\sigma^2+\theta^2\nu)t$. The parameter $\sigma$ alone is not total return volatility when $\theta\neq0$.

Why is omega needed?

$X_t$ has a non-zero exponential moment. The correction $\omega$ offsets it so that the discounted total-return stock is a martingale under the selected pricing measure.

Does VG have a unique risk-neutral measure?

No. The market is incomplete. A measure-change rule or direct option calibration is needed before applying the martingale correction.

How does VG differ from Merton jump diffusion?

Merton combines diffusion with a finite number of Poisson jumps. VG is pure jump with infinitely many small jumps on every interval.

Which papers introduced variance gamma?

Madan and Seneta's 1990 paper introduced the share-return model. Madan, Carr and Chang's 1998 paper developed the three-parameter option-pricing framework. The characteristic function is compact. The measure choice and term structure are not.