The Variance-Gamma Process: Returns With Skew and Kurtosis

Geometric Brownian motion produces lognormal returns. The lognormal is fine for many applications, but it has two well-known failures for financial data: symmetric tails (no skew) and insufficient kurtosis (not fat enough tails). Empirically, equity returns are negatively skewed (large losses more common than large gains) and have kurtosis well above 3 (the Gaussian value).

The variance-gamma (VG) process addresses both issues. It is a Lévy process with three parameters — controlling variance, kurtosis, and skew — and it has a tractable characteristic function. The VG model is one of the cleanest extensions of GBM that captures the empirical features of equity returns. This article derives its SDE, computes its moments, and shows why it fits the data better.

Why a New Process?

The empirical stylized facts of equity returns:

• Negative skew. Large losses occur more frequently than large gains. The left tail is heavier.
• High kurtosis. The distribution has fatter tails than Gaussian. Daily returns often have kurtosis of 5-15, well above 3.
• Volatility clustering. Volatility comes in waves. (Not addressed by static distributions.)
• Aggregational Gaussianity. Over long horizons, returns look more Gaussian. (VG also exhibits this.)

GBM produces lognormal returns, which have positive skew and kurtosis of $3 + (e^{\sigma^2 T} - 1)(e^{\sigma^2 T} + 2)$, which for typical $\sigma$ and $T$ is around 3.5-5 for daily data. Not bad, but skew is wrong (positive when empirical is negative), and tails are still too thin for short-horizon data.

The VG process produces returns that are negatively skewed and have arbitrary kurtosis (parameterized by $\nu$). The model fits daily and weekly equity returns much better than GBM.

The Variance-Gamma Process as a Time-Changed Brownian Motion

The VG process is a Brownian motion with stochastic time change. Specifically, let $W_t$ be a Brownian motion and $\gamma_t$ be a Gamma process with unit mean rate. Define

where $\theta$ and $\sigma$ are constants. This is the VG process.

The intuition: the Gamma process $\gamma_t$ is a positive, increasing Lévy process with Gamma-distributed increments. When $\gamma_t$ jumps, $X_t$ jumps. When $\gamma_t$ drifts slowly, $X_t$ drifts slowly. The result is a continuous-time Lévy process with skewed, fat-tailed increments.

The drift $\theta$ of the Gamma-time-changed Brownian motion controls skew. Positive $\theta$ gives positive skew (right tail heavy); negative $\theta$ gives negative skew (left tail heavy). For equity returns, we want $\theta < 0$.

The variance $\sigma^2$ controls the overall scale of the fluctuations. The variance of the Gamma process controls the variance of $X_t$. Three parameters: $\sigma$, $\theta$, and the variance $\nu$ of the Gamma process (sometimes written as $\nu = \text{Var}(\gamma_1)$).

The Three Parameters

The standard VG parameterization has parameters $(\sigma, \nu, \theta)$:

• $\sigma > 0$ is the volatility of the Brownian motion component.
• $\nu > 0$ is the variance of the Gamma time change.
• $\theta \in \mathbb{R}$ is the drift of the Brownian motion (and the skew parameter).

The marginal distribution of $X_1$ is the variance-gamma distribution, with density

where $K$ is the modified Bessel function of the second kind. This is a complicated closed form. In practice, you use the characteristic function for computation.

The mean of $X_1$ is $\theta$ (per unit of expected Gamma time). The variance is $\sigma^2 + \theta^2 \nu$. The skewness is

and the kurtosis is

For equity returns, $\theta < 0$ gives negative skew, and $\nu > 0$ gives excess kurtosis. The two together capture the empirical features.

The SDE Representation

The VG process has the SDE representation

with $\gamma_t$ a Gamma process. Equivalently, conditional on $\gamma_t = g$, $X_t \sim \mathcal{N}(\theta g, \sigma^2 g)$, i.e., $X_t$ is a mixture of Gaussians with Gamma-distributed variance.

This is the cleanest way to think about VG. It is a normal distribution with random variance, where the variance follows a Gamma distribution. The mixture of Gaussians with Gamma mixing produces a distribution with heavier tails than any single Gaussian.

For simulation, you can either:

1. Draw $\gamma_t$ from the Gamma distribution, then draw $X_t$ from $\mathcal{N}(\theta \gamma_t, \sigma^2 \gamma_t)$.
2. Use the series representation: $X_t = \theta \gamma_t + \sigma \sqrt{\gamma_t} Z$ where $Z \sim \mathcal{N}(0, 1)$ independent of $\gamma_t$.

Both work. The first is conceptually cleaner.

Characteristic Function

The characteristic function of the VG distribution is

This is one of the cleanest characteristic functions in finance. It exists in closed form, decays like a polynomial (heavy tails), and is easy to invert numerically.

For option pricing, you compute the characteristic function of $\log S_T$ and use Fourier inversion. The Carr-Madan formula or the COS method gives the option price in milliseconds.

The exponent $-t/\nu$ reflects the scaling: at time $t$, the characteristic function is the $t/\nu$-th power of the unit-time characteristic function. This is the Lévy property: VG is a Lévy process with stationary independent increments.

VG as a Model for Stock Prices

To model the stock price, take $S_t = S_0 \exp(X_t + \omega t)$, where $\omega$ is an adjustment to make the risk-neutral mean of $S_t$ equal to $S_0 e^{rt}$. The adjustment is

which makes the discounted stock $e^{-rt} S_t$ a martingale.

Under this model, the log-returns $\log(S_T / S_0)$ are VG-distributed (up to the linear drift $\omega T$). The distribution has the empirical features: negative skew and high kurtosis.

Option Pricing via Fourier Methods

The option price under VG is

Substituting $S_T = S_0 e^{X_T + \omega T}$ and using the characteristic function of $X_T$, we get

Using Fourier inversion,

where $\phi_{X_T}$ is the characteristic function of $X_T$ shifted. The integrals can be computed numerically with FFT-based methods.

The Carr-Madan (1999) approach computes the call price as the inverse Fourier transform of a damped characteristic function. The COS method (Fang and Oosterlee, 2008) uses a cosine expansion of the density. Both give prices accurate to 5-6 digits in milliseconds.

Calibration to Market Data

Calibrate $(\sigma, \nu, \theta)$ to fit the option prices observed in the market. The standard approach:

1. For a given parameter set, compute the option prices via Fourier methods.
2. Minimize the sum of squared differences between model and market prices.
3. The result is the calibrated VG parameter set.

The VG model produces an implied volatility surface with skew and smile, unlike Black-Scholes which gives flat vol. The skewness parameter $\theta$ controls the skew of the surface; $\nu$ controls the smile (curvature).

For a single maturity, VG has 3 parameters. The option surface at one maturity typically has hundreds of quotes. The fit is highly over-determined, and the calibration is stable.

import numpy as np
from scipy.optimize import minimize

def vg_char_func(u, t, sigma, nu, theta):
    """Characteristic function of VG at time t."""
    return (1 - 1j * u * theta * nu + 0.5 * sigma**2 * nu * u**2) ** (-t / nu)

def vg_call_price(S0, K, r, T, sigma, nu, theta):
    """VG call price via Fourier inversion (sketch)."""
    # Apply Carr-Madan formula with damping
    alpha = 1.5
    u = np.linspace(1e-5, 200, 10000)
    du = u[1] - u[0]
    log_K = np.log(K)
    cf = vg_char_func(u - 1j * (alpha + 1), T, sigma, nu, theta)
    # ... (omitted for brevity)
    return call_price

This is the skeleton of a VG pricer. The full implementation uses FFT or numerical integration of the characteristic function. The output is the call price as a function of $(\sigma, \nu, \theta)$. Calibration is then a least-squares optimization.

Comparison to Other Models

How does VG compare to other fat-tailed models?

VG vs. Merton jump-diffusion. Merton adds log-normal jumps to GBM. The result has fatter tails and can have skew. VG is similar in spirit but uses continuous-time stochastic time change rather than discrete jumps. VG is more analytically tractable (clean characteristic function); Merton has a sum-of-Black-Scholes formula.

VG vs. CGMY. CGMY (named after Carr, Geman, Madan, Yor) is a four-parameter Lévy process that generalizes VG. VG is the special case where the CGMY parameters satisfy $C = G = M = 0$ and $Y = 0$, leaving only three effective parameters. CGMY allows infinite activity of jumps (any time interval has infinitely many small jumps), while VG has finite activity per unit time.

VG vs. Heston stochastic vol. Heston models stochastic variance but assumes the underlying is Gaussian conditional on the variance. The result is a stochastic-vol model with closed-form characteristic function. VG is a static-distribution model (no stochastic vol), but with non-Gaussian marginals.

VG vs. SABR. SABR is a stochastic-vol model specifically designed to fit the volatility surface. It has closed-form approximation (Hagan’s formula) for the implied vol. VG fits the same data via the Fourier approach.

For fitting the volatility surface with a clean, fast model, VG is a strong choice. For capturing path-dependent features (volatility clustering, leverage effect), you need stochastic vol or regime switching.

Why VG Captures Equity Returns Well

The empirical fit of VG to equity returns is excellent. The Madan-Seneta (1990) and Madan-Milne (1991) papers showed that VG fits daily S&P 500 returns much better than GBM. The negative skew is captured by $\theta < 0$. The kurtosis is captured by $\nu$. The volatility is captured by $\sigma$.

The key feature: VG is a mixture of Gaussians with Gamma mixing. This is a flexible distributional family. Many empirical distributions can be approximated by such mixtures. VG is one of the simplest mixtures with closed-form characteristic function.

For longer-horizon returns, VG also exhibits aggregational Gaussianity. The VG distribution at time $T$ has kurtosis $3 + 3/T$ (in the symmetric case with $\theta = 0$), which approaches 3 (Gaussian) as $T \to \infty$. This matches the empirical observation that long-horizon returns look more normal than short-horizon returns.

Simulation of VG Paths

For Monte Carlo pricing of path-dependent options, simulate VG paths. The simplest approach:

1. Draw Gamma-distributed time increments $\Delta \gamma_i \sim \text{Gamma}(\Delta t / \nu, \nu)$.
2. Draw normal increments $Z_i \sim \mathcal{N}(0, 1)$ independent.
3. Set $\Delta X_i = \theta \Delta \gamma_i + \sigma \sqrt{\Delta \gamma_i} Z_i$.
4. Cumulate: $X_t = \sum_{i} \Delta X_i$.

The stock price is $S_t = S_0 \exp(X_t + \omega t)$. This gives a path with VG-distributed increments and the right martingale property.

For path-dependent options (Asian, lookback, barrier), the full path matters. The simulation captures the path correctly because the increments are i.i.d. VG-distributed. Just be careful to use enough time steps to capture the path features.

import numpy as np

def simulate_vg_paths(S0, T, n_steps, n_paths, sigma, nu, theta, r, seed=42):
    rng = np.random.default_rng(seed)
    dt = T / n_steps
    omega = np.log(vg_martingale_correction(T, sigma, nu, theta, r)) / T
    # Draw Gamma increments
    shape = dt / nu
    scale = nu
    d_gamma = rng.gamma(shape, scale, size=(n_paths, n_steps))
    # Draw Normal increments
    Z = rng.standard_normal((n_paths, n_steps))
    # VG increments
    dX = theta * d_gamma + sigma * np.sqrt(d_gamma) * Z
    # Cumulative sum
    X = np.cumsum(dX, axis=1)
    t = np.linspace(dt, T, n_steps)
    S = S0 * np.exp(X + omega * t)
    return np.concatenate([np.full((n_paths, 1), S0), S], axis=1)

This produces paths of $S_t$ with VG-distributed log-returns. The martingale correction $\omega$ ensures risk-neutral pricing.

The Madan-Seneta Result

Madan and Seneta (1990) showed that the VG process has an elegant connection to the symmetric hyperbolic distribution. The hyperbolic distribution is a four-parameter family with closed-form density. VG is a special case.

The hyperbolic family includes:

• Normal ($\nu \to 0$).
• VG (finite $\nu > 0$).
• Hyperbolic (full four parameters).
• Generalized hyperbolic (further extensions).

For practical work, VG is the most useful because it has only three parameters and a clean characteristic function. The extra parameters of the hyperbolic family are often not needed to fit equity return data.

A Note on Numerical Stability

The VG characteristic function involves powers of a complex number. For large $|u|$, this can overflow or underflow. Implementations should use log-space arithmetic:

The log of a complex number requires care with branch cuts. Most implementations use numpy.log of complex numbers and handle the imaginary part carefully.

For option pricing, the integrals in the Fourier inversion can be sensitive to the choice of integration contour. The Carr-Madan damping factor $\alpha$ must be chosen to ensure the integrals converge. Typical values are $\alpha = 1.0$ to $2.0$.

What to Take Away

The variance-gamma process is a three-parameter Lévy process: $\sigma$ controls volatility, $\nu$ controls kurtosis, $\theta$ controls skew. It is a Brownian motion with Gamma-distributed stochastic time change. The characteristic function is closed-form, the density is in terms of Bessel functions, and option prices can be computed by Fourier methods in milliseconds.

For equity returns, VG captures the negative skew and excess kurtosis that GBM misses. The model is widely used for pricing equity options, FX options, and some credit derivatives. It is the cleanest fat-tailed, skewed extension of GBM with a tractable characteristic function.

References: Madan and Seneta (1990), Journal of Business; Madan and Milne (1991); Carr, Geman, Madan, and Yor (2002); Cont and Tankov (2004), Financial Modelling with Jump Processes, Chapter 4. The VG paper has over 1,500 citations and is one of the standard references for Lévy process modeling in finance. Read Cont and Tankov for the broader Lévy framework. We cover Lévy processes in general in the next article.

Why Bessel Functions Show Up

The density of the VG distribution involves the modified Bessel function of the second kind, $K_\alpha(z)$. This function appears throughout applied mathematics whenever you have mixtures of Gaussians with Gamma-distributed variances. The mathematical reason: the integral

is a standard Laplace-type identity, and the VG density comes from exactly this integral with $g$ being the Gamma mixing variable.

For numerical evaluation, Bessel functions are built into scipy.special and most numerical libraries. They are well-behaved for moderate arguments and require care for very large or very small arguments (use asymptotic expansions).

For option pricing, you typically do not need to evaluate the density directly. The characteristic function is enough. But if you need the density for risk management or stress testing, the Bessel form is the way to go.

VG and the Self-Decomposable Property

The VG process is self-decomposable: for any $0 < s < t$, $X_t$ has the same distribution as $X_s +$ an independent increment with VG distribution. This is the Lévy property — independent, stationary increments. It is what makes the VG process a Lévy process.

Self-decomposability is a stronger property than just having stationary independent increments. It says that the distribution at any time can be written as a sum of an earlier distribution plus an independent increment. VG satisfies this. Many other Lévy processes do too.

For pricing, self-decomposability is useful because the conditional characteristic function $\mathbb{E}[e^{iu(X_t - X_s)} | X_s] = \phi_X(u)^{(t-s)/t \cdot 1}$ factorizes. This makes Fourier methods clean and efficient.

VG and Infinite Divisibility

VG is also infinitely divisible: for any integer $n$, the distribution of $X_t$ can be written as the sum of $n$ i.i.d — random variables. Specifically, $X_t \overset{d}{=} \sum_{i=1}^n X_{t/n}$. This is the Lévy-Khintchine representation in another guise.

For Monte Carlo, infinite divisibility lets you simulate VG increments over arbitrarily small time intervals by breaking the increment into many small pieces. This is useful for path-dependent options where you need fine time resolution.

The Gamma Process

The Gamma process $\gamma_t$ that drives the VG time change has Gamma-distributed increments: $\gamma_{t+s} - \gamma_t \sim \text{Gamma}(s/\nu, \nu)$, with shape $s/\nu$ and scale $\nu$. The mean is $s$ and the variance is $s\nu$. For small $s$, the increments are small and Gamma-like (with low variance). For large $s$, the increments are larger and more variable.

The Gamma process has paths of bounded variation. Each increment is positive. So $\gamma_t$ is strictly increasing. This means $W_{\gamma_t}$ runs through Brownian motion at a variable speed: faster when $\gamma_t$ jumps up, slower when it does not.

For simulation, draw $\Delta \gamma_i \sim \text{Gamma}(\Delta t / \nu, \nu)$ for each time step. The increments are independent. Cumulate to get $\gamma_t$.

VG and the Laplace Distribution

In the limit $\nu \to 0$, the VG distribution becomes Gaussian. In the limit $\nu \to \infty$, the VG distribution becomes the Laplace distribution (double exponential), which is a heavy-tailed distribution with exponential decay. The Laplace has fatter tails than Gaussian but thinner than stable distributions.

The parameter $\nu$ interpolates between these two limits. Small $\nu$ gives near-Gaussian behavior. Large $\nu$ gives heavy-tailed behavior. For equity returns, the calibrated $\nu$ is typically around 0.1-0.5 (per year), giving moderate tails.

The Laplace distribution has the density $f(x) = (1/2b) e^{-|x - \mu|/b}$, which is exponential in $|x - \mu|$. The tail probability decays like $e^{-|x|/b}$, which is faster than polynomial but slower than Gaussian’s $e^{-x^2/(2\sigma^2)}$. The kurtosis of the Laplace is 6.

Pricing Exotic Options under VG

For path-dependent options (Asian, lookback, barrier), VG requires simulation rather than closed-form Fourier methods. The standard approach:

1. Simulate VG paths using the Gamma-time-changed Brownian motion approach.
2. Compute the payoff on each path.
3. Average discounted payoffs for the price.

The simulation is exact because VG is a Lévy process with i.i.d — increments. Just draw the right number of time steps and accumulate.

For variance reduction, use antithetic variates (draw $Z$ and $-Z$ for paired paths) or control variates (use the GBM price as a control). Both work well for VG because the model is close enough to GBM that the control is highly correlated.

For barrier options, the simulation requires care: the process jumps (because of the Gamma time change), so a barrier that is touched by a jump may be missed by a fine-time-step simulation. Use small enough time steps or analytic corrections for the jump probability.

A Comparison: VG Skew vs. Vol Smile

The VG model produces an implied volatility skew. The skew is downward-sloping (negative) for $\theta < 0$, which is what equity markets show. The slope of the skew at-the-money is approximately

Higher $|\theta|$ gives a steeper skew. Higher $\sigma$ gives a flatter skew. Longer maturity $T$ gives a flatter skew (the skew decays as $1/\sqrt{T}$).

This is consistent with empirical observations: short-dated options have steeper skews than long-dated options. VG captures this naturally. The Heston model also captures it but via stochastic vol rather than static skew.

For risk management, the VG skew parameter $\theta$ is approximately constant across time (unlike the Heston parameters which may drift). This makes VG easier to use for stress testing and scenario analysis.

VG and the Moment Explosion Problem

VG has finite moments of all orders. The tails decay fast enough that all powers of $X$ are integrable. This is in contrast to stable distributions (covered in the next article), which have infinite moments beyond a critical order.

For option pricing, finite moments mean that the call price has a finite expectation. For risk management, finite moments mean that Value-at-Risk and Expected Shortfall are well-defined. These are practical advantages of VG over heavier-tailed alternatives.

For very deep out-of-the-money options or stress scenarios, the VG model may underestimate tail risk. For those applications, you need heavier-tailed models (CGMY, stable, or jump-diffusion with power-law jumps). For typical option pricing, VG is sufficient.

The VG Model in Industry

The VG model is used by several major banks for equity and FX option pricing. It is particularly popular for pricing options on assets that have heavy tails but are not subject to extreme jumps (like single stocks in stable markets). The clean characteristic function and fast Fourier pricing make it efficient for production use.

For multi-asset options (basket options, options on the minimum of a set of stocks), VG has been extended to multi-asset versions with correlated Gamma processes. The characteristic function is still tractable, though more complex.

The main competitor is the CGMY model, which has four parameters and can capture even heavier tails. CGMY is more flexible but harder to calibrate. VG is the sweet spot of tractability and flexibility for many applications.

An Alternative: The Normal Inverse Gaussian

The Normal Inverse Gaussian (NIG) distribution is another three-parameter Lévy process that fits equity returns well. It is the special case of the generalized hyperbolic distribution with $\lambda = -1/2$. The NIG has semi-heavy tails (heavier than Gaussian, lighter than stable) and analytical tractability.

VG and NIG are both popular for option pricing. VG has cleaner moment formulas; NIG has a slightly simpler density. For most applications, the choice between them comes down to which fits your data better. Both are reasonable choices for fat-tailed equity return modeling.

The Barndorff-Nielsen (1997) paper introduced NIG. The original NIG paper is short and well-written. Read it if you want to compare VG and NIG side-by-side.