Stochastic Volatility

[math] \newcommand{\indicator}[1]{\mathbbm{1}_{\left[ {#1} \right] }} \newcommand{\Real}{\hbox{Re}} \newcommand{\HRule}{\rule{\linewidth}{0.5mm}} \newcommand{\mathds}{\mathbb}[/math]

Stochastic volatility is another popular way to fit the implied volatility smile. It has the advantage of assuming volatility brings another source of randomness, and this randomness can be hedge if it is also a Brownian motion. Furthermore, calibration of stochastic volatility models does have some robustness to the market's daily changes. However, whereas local volatility models are generally calibrated to fit the implied volatility surface (i.e. to fit options of all maturities), stochastic volatility models have trouble fitting more than one maturity at a time. Consider the following stochastic volatility model,

[[math]] \begin{eqnarray} \label{eq:SVM} dS_t&=&\mu S_tdt+\sigma(X_t)S_tdW_t\\ \label{eq:dX} dX_t &=&\alpha(X_t)dt+\beta(X_t)dB_t \end{eqnarray} [[/math]]

where [math]W_t[/math] and [math]B_t[/math] are Brownian motions with correlation [math]\rho\in(-1,1)[/math] such that [math]\mathbb E[dW_tdB_t]=\rho dt[/math], [math]\sigma(x) \gt 0[/math] is the volatility function, and [math]X_t[/math] is the volatility process and is usually mean-reverting. An example of a mean-revering process is the square-root process,

[[math]]dX_t = \kappa(\bar X-X_t)dt+\gamma\sqrt{X_t}dB_t\ .[[/math]]

The square-root process is used in the Heston model (along with [math]\sigma(x) =\sqrt{x}[/math]), and usually relies on what is known as the Feller condition ([math]\gamma^2\leq 2\bar X\kappa[/math]) so that [math]X_t[/math] never touches zero.

Hedging and Pricing

Assuming that [math]\alpha[/math] and [math]\beta[/math] are ‘well-behaved’ functions, we define the differential operator

[[math]]\mathcal L \doteq\frac{\partial}{\partial t}+\frac{\sigma^2(x)s^2}{2}\frac{\partial^2}{\partial s^2}+\mu s \frac{\partial}{\partial s}+\frac{\beta^2(x)}{2}\frac{\partial^2}{\partial x^2}+\alpha(x)\frac{\partial}{\partial x}+\rho\beta(x)\sigma(x)s\frac{\partial^2}{\partial x\partial s}\ ; [[/math]]

the interpretation of [math]\mathcal L[/math] is the following: [math]\frac{\sigma^2(X_t)S_t^2}{2}\frac{\partial^2}{\partial s^2}+\mu S_t \frac{\partial}{\partial s}[/math] are the generator of [math]S_t[/math]; [math]\frac{\beta^2(X_t)}{2}\frac{\partial^2}{\partial x^2}+\alpha(X_t)\frac{\partial}{\partial x}[/math] are the generator of [math]X_t[/math]; [math]\rho\beta(X_t)\sigma(X_t)S_t\frac{\partial^2}{\partial x\partial s}[/math] is the cross-term; and the remaining are the separate stochastic terms. We use when applying bi-variate It\^o Lemma for the process in \eqref{eq:SVM} and \eqref{eq:dX}:\\ For any twice differentiable function [math]f(t,s,x)[/math], the It\^o differential is

[[math]] \begin{align} \label{eq:2Dito} df(t,S_t,X_t) &= \mathcal Lf(t,S_t,X_t)dt+\sigma(X_t)S_t\frac{\partial}{\partial s}f(t,S_t,X_t)dW_t+\beta(X_t)\frac{\partial}{\partial x}f(t,S_t,X_t)dB_t\ . \end{align} [[/math]]

The price [math]C(t,s,x)[/math] is the price of claim [math]\psi(s,x)[/math] paid at time [math]T[/math] given [math]S_t=s[/math] and [math]X_t=x[/math]. We hedge with a self-financing portfolio [math]V_t[/math] consisting of risky-asset, the risk-free bank account, and a 2nd derivative [math]C'[/math] that is settled at time [math]T' \gt T[/math]. This method for hedging stochastic volatility is called the Hull-White method. The self-financing condition is

[[math]]dV_t = a_tdS_t+b_tdC'(t,S_t,X_t)+r(V_t - a_tS_t-b_tC'(t,S_t,X_t))dt\ .[[/math]]

Now, taking a long position in [math]V_t[/math] and a short position in [math]C(t,S_t,X_t)[/math], we apply the bivariate It\^o lemma of \eqref{eq:2Dito} to get the dynamics of this new portfolio:

[[math]]d(V_t-C(t,S_t,X_t))[[/math]]

[[math]]=a_t(\mu S_tdt+\sigma(X_t)S_tdW_t)[[/math]]

[[math]]+ b_t\mathcal LC'(t,S_t,X_t)dt+b_t\sigma(X_t)S_t\frac{\partial}{\partial s}C'(t,S_t,X_t)dW_t+b_t\beta(X_t)\frac{\partial}{\partial x}C'(t,S_t,X_t)dB_t[[/math]]

[[math]]+r(V_t - a_tS_t-b_tC'(t,S_t,X_t))dt[[/math]]

[[math]]-\mathcal LC(t,S_t,X_t)dt-\sigma(X_t)S_t\frac{\partial}{\partial s}C(t,S_t,X_t)dW_t-\beta(X_t)\frac{\partial}{\partial x}C(t,S_t,X_t)dB_t\ .[[/math]]

We now choose [math]a_t[/math] and [math]b_t[/math] so that the [math]dW_t[/math] and the [math]dB_t[/math] terms vanish, which means that

[[math]] \begin{align*} &\sigma(X_t)S_t\frac{\partial}{\partial s}C(t,S_t,X_t) = a_t\sigma(X_t)S_t+ b_t\sigma(X_t)S_t\frac{\partial}{\partial s}C'(t,S_t,X_t)\\ &\\ &\beta(X_t)\frac{\partial}{\partial x}C(t,S_t,X_t)=b_t\beta(X_t)\frac{\partial}{\partial x}C'(t,S_t,X_t) \end{align*} [[/math]]

or

[[math]] \begin{align} \label{eq:a_t} &b_t = \frac{\frac{\partial}{\partial x}C(t,S_t,X_t)}{\frac{\partial}{\partial x}C'(t,S_t,X_t)}\\ \label{eq:b_t} &a_t = \frac{\partial}{\partial s}C(t,S_t,X_t)-b_t\frac{\partial}{\partial s}C'(t,S_t,X_t)\ . \end{align} [[/math]]

Plugging these solutions into the differential, the [math]\alpha(X_t)\frac{\partial}{\partial x}[/math] terms cancel, the [math]\mu S_t\frac{\partial}{\partial s}[/math] terms cancel, and by the same arbitrage argument as Black-Scholes, it follows that [math]V_t-C(t,S_t,X_t)[/math] must grow at the risk-free rate:

[[math]] \begin{align*} &d(V_t-C_t(t,S_t,X_t))\\ &= a_t\mu S_t dt+b_t\mathcal LC'(t,S_t,X_t)dt+r(V_t - a_tS_t-b_tC'(t,S_t,X_t))dt -\mathcal LC(t,S_t,X_t)dt\\ &=r(V_t-C(t,S_t,X_t))dt\ . \end{align*} [[/math]]

Inserting [math]a_t[/math] and [math]b_t[/math] from equations \eqref{eq:a_t} and \eqref{eq:b_t}, and after some rearranging of terms (and dropping the [math]dt[/math]'s), we get

[[math]]\left(\mathcal L-(\mu-r)S_t\frac{\partial}{\partial s}-r\right)C(t,S_t,X_t)= b_t\left(\mathcal L-(\mu-r)S_t\frac{\partial}{\partial s}-r\right)C'(t,S_t,X_t)[[/math]]

and dividing both sides by [math]\frac{\partial}{\partial x}C(t,S_t,X_t)[/math], we get

[[math]]\frac{\left(\mathcal L-(\mu-r)S_t\frac{\partial}{\partial s}-r\right)C(t,S_t,X_t)}{\frac{\partial}{\partial x}C(t,S_t,X_t)}= \frac{\left(\mathcal L-(\mu-r)S_t\frac{\partial}{\partial s}-r\right)C'(t,S_t,X_t)}{\frac{\partial}{\partial x}C'(t,S_t,X_t)}\doteq R(t,s,x)\ .[[/math]]

The left-hand side of this equation does not depend on specifics of [math]C'[/math] (e.g. maturity and strikes), and the right-hand side does not depend on specifics of [math]C[/math]. Therefore, there is a function [math]R(t,s,x)[/math] that does not depend on parameters such as strike price and maturity, and it equates the two sides of the expression. We let [math]R[/math] take the form

[[math]]R(t,s,x) = -\alpha(x)+\beta(x)\Lambda(t,s,x)[[/math]]

where

[[math]]\Lambda(t,s,x) = \rho\frac{\mu-r}{\sigma(x)}+g(t,s,x)\sqrt{1-\rho^2}\ ,[[/math]]

where [math]g[/math] is an arbitrary function. Hence, we arrive at the PDE for the price [math]C(t,s,x)[/math] of a European Derivative with stochastic vol:

Proposition (Pricing PDE for Stochastic Volatility)

Given the stochastic volatility model of equations \eqref{eq:SVM} and \eqref{eq:dX}, the price of a European claim with payoff [math]\psi(S_T)[/math] satisfies

[[math]] \begin{align} \nonumber \Bigg(\frac{\partial}{\partial t}+\frac{\sigma^2(x)s^2}{2}\frac{\partial^2}{\partial s^2}+\rho\beta(x)\sigma(x)s\frac{\partial^2}{\partial x\partial s}+\alpha(x)\frac{\partial}{\partial x}+\frac{\beta^2(x)}{2}\frac{\partial^2}{\partial x^2}+rs\frac{\partial}{\partial s}-r\Bigg)C&\\ \label{eq:SVpde} =\beta(x)\Lambda(t,s,x)\frac{\partial}{\partial x}C& \end{align} [[/math]]

with [math]C\Big|_T=\psi(s,x)[/math].

The PDE in \eqref{eq:SVpde} can be solved with Feynman-Kac,

[[math]]C(t,s,x) = e^{-r(T-t)}\mathbb E^Q\{\psi(S_T,X_T)|S_t=s,X_t=x\}\ .[[/math]]

which yields a description of an EMM, [math]\mathbb Q[/math], under which the stochastic volatility model is

[[math]] \begin{eqnarray*} \frac{dS_t}{S_t}&=&r dt+\sigma(X_t)dW_t^Q\\ dX_t&=&\left(\alpha(X_t)-\beta(X_t)\Lambda(t,S_t,X_t)\right)dt+\beta(X_t)dB_t^Q \end{eqnarray*} [[/math]]

where [math]W_t^Q[/math] and [math]B_t^Q[/math] are [math]\mathbb Q[/math]-Brownian motions with the same correlation [math]\rho[/math] as before. The interpretation of [math]\Lambda(t,s,x)[/math] is that it is the market price of volatility risk. One can think of a market as complete when there are as many non-redundant assets as there are sources of randomness. In this sense we have ‘completed’ this market with stochastic volatility by including the derivative [math]C'[/math] among the traded assets. This has allowed us to replicate the European claim and obtain a unique no-arbitrage price that must be the solution to \eqref{eq:SVpde}.\\

Example

[[math]] \begin{align} \label{eq:dSheston} &dS_t=r S_tdt+\sqrt{X_t}dW_t^Q\\ \label{eq:dXheston} &dX_t=\kappa(\bar X-X_t)dt+\lambda X_tdt+\gamma\sqrt{X_t}dB_t^Q \end{align} [[/math]]

where [math]\lambda X_t[/math] is the assumed form of the volatility price of risk. This model can be re-written,

[[math]]dX_t = \tilde\kappa(\tilde{\bar X}-X_t)dt+\gamma\sqrt{X_t}dB_t^Q[[/math]]

where [math]\tilde\kappa =\kappa-\lambda[/math] and [math]\tilde{\bar X} = \frac{\kappa\bar X}{\kappa-\lambda}[/math]. It is important to have the Feller condition

[[math]]\gamma^2\leq 2\tilde\kappa\tilde{\bar X}[[/math]]

otherwise the PDE requires boundaries for the event [math]X_t = 0[/math]. The process [math]X_t[/math] is degenerate, but it is well-known that the Feller condition is the critical assumption for the application of Feynman-Kac. The fit of the Heston model to the implied volatility of market data is shown in Figure.

The fit of the Heston model to implied volatility of market data on July 27 of 2012, with different maturities. Notice that the Heston model doesn't fit very well to the options with shortest time to maturity.

Monte Carlo Methods

A slow but often reliable method for pricing derivatives is to approximate the risk-neutral expectation with independent samples. Essentially, sample paths of Brownian motion are obtained using a pseudo-random number generated, the realized derivative payoff is computed for each sample, and then the average of these sample payoffs is an approximation of the true average. For example, let [math]\ell=1,2,\dots,N[/math], let [math]S_T^{(\ell)}[/math] be the index for an independent sample from [math]S_T[/math]'s probability distribution. If [math]\mathbb E^QS_T^2 \lt \infty[/math], then by the law of large numbers we have

[[math]]\frac 1N\sum_{\ell=1}^N(S_T^{(\ell)}-K)^+\stackrel{N\rightarrow\infty}{\longrightarrow}\mathbb E^Q(S_T-K)^+\ .[[/math]]

For stochastic differential equations, the typical method is to generate samples that are approximately from the same distribution as the price. The basic idea is to use an Euler-Maruyama scheme,

[[math]] \begin{equation} \label{eq:EMscheme} \log \tilde S_{t+\Delta t}^{(\ell)} = \log \tilde S_t^{(\ell)}+\left(\mu-\frac 12\sigma^2\right)\Delta t+\sigma(W_{t+\Delta t}^{(\ell)}-W_t^{(\ell)})\ , \end{equation} [[/math]]

for some small [math]\Delta t \gt 0[/math], with

[[math]]W_{t+\Delta t}^{(\ell)}-W_t^{(\ell)}\sim iid \mathcal N(0,\Delta t)\ .[[/math]]

This scheme is basically the discrete backward Riemann sum from which the It\^o integral was derived (see Chapter). In the limit, the Monte-Carlo average will be with little-o of the true expectation,

[[math]]\lim_{N\rightarrow\infty} \frac 1N\sum_{\ell=1}^N(\tilde S_T^{(\ell)}-K)^+ = \mathbb E^Q(S_T-K)^++o(\Delta t)\ .[[/math]]

The scheme in \eqref{eq:EMscheme} can easily be adapted for a local volatility model, and can also be adapted to stochastic volatility models, but the latter will also need a scheme for generating the volatility process [math]X_t[/math]. Example

[[math]] \begin{align*} &\frac{dS_t}{S_t} = rdt+e^{X_t}dW_t^Q\\ &dX_t = \kappa(\bar X-X_t)dt+\gamma dB_t^Q \end{align*} [[/math]]

with correlation [math]\rho[/math] between [math]W_t^Q[/math] and [math]B_t^Q[/math]. The Euler-Maruyama scheme is

[[math]] \begin{align*} &\log \tilde S_{t+\Delta t} = \log \tilde S_t+\left(r-\frac 12e^{2\tilde X_t}\right)\Delta t+e^{\tilde X_t}\left(\rho (B_{t+\Delta t}-B_t)+\sqrt{1-\rho^2}(W_{t+\Delta t}-W_t)\right)\\ &\tilde X_{t+\Delta t} = \tilde X_t+\kappa\left(\bar X-\tilde X_t \right)\Delta t+\gamma(B_{t+\Delta t}-B_t) \end{align*} [[/math]]

where [math]B_{t+\Delta t}^{(\ell)}-B_t^{(\ell)}[/math] and [math]W_{t+\Delta t}^{(\ell)}-W_t^{(\ell)}[/math] are increments of independent Brownian motions.

Sometimes, there more clever methods of sampling need to be employed. For instance, the Heston model: Example

[[math]]\tilde X_{t+\Delta t}=\tilde X_t+\kappa(\bar X-\tilde X_t)\Delta t+\gamma\sqrt{\tilde X_t}\left(B_{t+\Delta t}-B_t\right)[[/math]]

will allow the volatility process to become negative. An alternative is to consider the Stratonovich-type differential,

[[math]]dX_t = \left(\kappa(\bar X-X_t)-\frac{\gamma^2}{2}\right)dt+\gamma\sqrt{X_t}\circ dB_t[[/math]]

and then use an implicit Euler-Maruyama scheme

[[math]]\tilde X_{t+\Delta t}-\tilde X_t =\left( \kappa(\bar X-\tilde X_{t+\Delta t})-\frac{\gamma^2}{2}\right)\Delta t+\gamma \sqrt{\tilde X_{t+\Delta t}} \left(B_{t+\Delta t}-B_t\right)\ .[[/math]]

which can be solved for [math]\sqrt{\tilde X_{t+\Delta t}}[/math] using the quadratic equation. Notice that a condition for the discriminant to be real is the Feller condition that [math]\gamma^2\leq 2\kappa\bar X[/math].

In general, Monte Carlo methods are good because they can be used to price pretty much any type of derivative instrument, so long as the EMM is given. For instance, it is very easy to price an Asian option with Monte Carlo. Another advantages of Monte Carlo is that there is no need for an explicit solution to any kind of PDE; the methodology for pricing is to simply simulate and average. However, as mentioned above, this can be very slow, and so there is a need (particularly among practitioners) for faster ways to compute prices. Another important topic to remark on is importance sampling, which refers to simulations generated under an equivalent measure for which certain rare events are more common. For instance, suppose one wants to use Monte Carlo to price a call option with extreme strike. It will be a waste of time to generate samples from the EMM's model if the overwhelming majority of them will finish out of the money. Instead, one can sample from another model that has more volatility, so that more samples finish in the money. Then, assign importance weights based on ratio of densities and take the average. See for further reading on Monte Carlo methods in finance.

Fourier Transform Methods

The most acclaimed way to compute a European option price is with the so-called Fourier transform methods. In particular, the class of affine models in have tractable pricing formulae because they have explicit expressions for their Fourier transforms. The essential ingredient to these methods is the characteristic function,

[[math]]\phi_T(u)\doteq\mathbb E^Q\exp(iu\log(S_T))\qquad\forall u\in\mathbb C,~i=\sqrt{-1}\ ,[[/math]]

where [math]\mathbb E^Q[/math] now denotes expectation under an EMM. The function [math]\phi_T(u)[/math] is the Fourier transform of the [math]\log S_T[/math]'s distribution function,

[[math]]\phi_T(u) = \int_{-\infty}^\infty e^{iu s}f(s)ds[[/math]]

where [math]f(s)[/math] is the density function for [math]\log S_T[/math]. The inverse Fourier transform is [math]f(s) = \tfrac{1}{2\pi}\int_{-\infty}^\infty e^{-iu s}\phi_T(u)du[/math], which leads to the cumulative distribution function

[[math]]F(s)-F(0)=\int_0^sf(v)dv = \frac{1}{2\pi}\int_{-\infty}^\infty \frac{1-e^{-ius}}{iu}\phi_T(u)du\ .[[/math]]

The crucial piece of information is the characteristic function, as any European claim can be priced provided there is also a Fourier transform for the payoff function. Let [math]p_T[/math] denote the density function on the risk-neutral distribution of [math]\log(S_T)[/math] and let [math]\widehat\psi[/math] denote the Fourier transform of the payoff function. It is shown in that

[[math]]\mathbb E^Q\psi(\log(S_T)) = \int_{-\infty}^\infty\psi(s)f(s)ds=\int_{-\infty}^\infty\frac{1}{2\pi} \int_{iz-\infty}^{iz+\infty}e^{-ius}\widehat \psi(u) du f(s)ds[[/math]]

[[math]]=\frac{1}{2\pi}\int_{iz-\infty}^{iz+\infty}\widehat \psi(u)\int_{-\infty}^\infty e^{-ius} p_T(s)ds du=\frac{1}{2\pi}\int_{iz-\infty}^{iz+\infty}\widehat \psi(u)\phi_T(-u) du[[/math]]

where the constant [math]z\in\mathbb R[/math] in the limit of integration is the appropriate strip of regularity for non-smooth payoff functions (see Table). This usage of characteristic functions is far reaching, as the class of L\'evy jump models are defined in terms of their characteristic function, or equivalently with the L\'evy-Khintchine representation (for more on financial models with jumps see ).

Fourier Transforms of Various Payoffs, [math]\widehat\psi(u)=\int_{-\infty}^\infty e^{ius}\psi(s)ds[/math] and [math]\psi(s)=\frac{1}{2\pi} \int_{iz-\infty}^{iz+\infty}e^{-ius}\widehat \psi(u) du[/math], [math]z=\Im(u)[/math].

asset	[math]\psi(s)[/math]	[math]\widehat\psi(u)[/math]	regularity strip
call	[math](e^s-K)^+[/math]	[math]-\frac{K^{iu+1}}{u^2-iu}[/math]	[math]z \gt 1[/math]
put	[math](K-e^s)^+[/math]	[math]-\frac{K^{iu+1}}{u^2-iu}[/math]	[math]z \lt 0[/math]
covered call; cash secured put	[math]\min(e^s,K)[/math]	[math]\frac{K^{iu+1}}{u^2-iu}[/math]	[math]0 \lt z \lt 1[/math]
cash or nothing call	[math]\indicator{e^s\geq K}[/math]	[math]-\frac{K^{iu}}{iu}[/math]	[math]z \gt 0[/math]
cash or nothing put	[math]\indicator{e^s\leq K}[/math]	[math]\frac{K^{iu}}{iu}[/math]	[math]z \lt 0[/math]
asset or nothing call	[math]e^s\indicator{e^s\geq K}[/math]	[math]-\frac{K^{iu+1}}{iu+1}[/math]	[math]z \gt 1[/math]
asset or nothing put	[math]e^s\indicator{e^s\leq K}[/math]	[math]\frac{K^{iu+1}}{iu+1}[/math]	[math]z \lt 0[/math]
Arrow-Debreu	[math]\delta(s-\log(K))[/math]	[math]K^{iu}[/math]	[math]z\in\mathbb R[/math]

Call-Option Pricing

An important tool is the formula of Gil-Pelaez:

Proposition (Gil-Pelaez Inversion Theorem)

The distribution function [math]F(s)[/math] is given by

[[math]] \begin{equation} \label{eq:gilPelaezInversion} F(s) = \frac12+\frac{1}{2\pi}\int_0^\infty \frac{e^{ius}\phi_T(-u)-e^{-ius}\phi_T(u)}{iu}du\ , \end{equation} [[/math]]

for all [math]s\in(-\infty,\infty)[/math].

Show Proof

First let's prove an identity,

[[math]]\mbox{sign}(v-s) = \frac2\pi\int_0^\infty \frac{\sin(u(v-s))}{u}du = \left\{\begin{array}{cc}-1&\mbox{if }v \lt s\\0&\mbox{if }v=s\\1&\mbox{if }v \gt s\end{array}\right. \ .[[/math]]

This integral is shown using residue calculus:

[[math]] \begin{align*} \int_0^\infty \frac{\sin(u)}{u}du &= \frac12\int_{-\infty}^\infty \frac{\sin(u)}{u}du=\frac12 \lim_{\epsilon\rightarrow 0}\lim_{R\rightarrow 0}\mbox{Im}\left(\int_{-R}^{-\epsilon}\frac{e^{iz}}{z}dz+\int_\epsilon^R\frac{e^{iz}}{z}dz\right)\ . \end{align*} [[/math]]

A closed path can be constructed from [math]-R[/math] to [math]R[/math], making a semi-circle of radius [math]\epsilon[/math] around zero, and then from [math]R[/math] back to [math]-R[/math] on a semi-circle of radius [math]R[/math]. Along this path the integral of [math]e^{iz}/z[/math] is zero,

[[math]]\int_{-R}^{-\epsilon}\frac{e^{iz}}{z}dz+\int_\epsilon^R\frac{e^{iz}}{z}dz=-\int_{C_\epsilon}\frac{e^{iz}}{z}dz-\int_{C_R}\frac{e^{iz}}{z}dz=0\ ,[[/math]]

where [math]C_\epsilon =\{\epsilon e^{i(\theta-\pi)}|0\leq \theta\leq \pi\}[/math] and [math]C_R=\{Re^{i\theta}|0\leq \theta\leq \pi\}[/math]. Now,

[[math]] \begin{align*} \left|\int_{C_R}\frac{e^{iz}}{z}dz \right|&=\frac1R\left|\int_0^\pi\frac{e^{iR\cos(\theta)-R\sin(\theta)}}{e^{i\theta}}d\theta\right|\\ &\leq \frac{1}{R}\int_0^\pi\left|\frac{e^{iR\cos(\theta)-R\sin(\theta)}}{e^{i\theta}}\right|d\theta \\ &= \frac{1}{R}\int_0^\pi e^{-R\sin(\theta)}d\theta \\ & \lt \frac{\pi}{R}\rightarrow 0\qquad\hbox{as }R\rightarrow\infty \ , \end{align*} [[/math]]

and

[[math]]\int_{C_\epsilon}\frac{e^{iz}}{z}dz = \int_\pi^0 \frac{1+O(\epsilon)}{\epsilon e^{i\theta}}i\epsilon e^{i\theta}d\theta = -i\pi +O(\epsilon)\ .[[/math]]

Hence,

[[math]]\int_0^\infty \frac{\sin(u)}{u}du = \frac12 \lim_{\epsilon\rightarrow 0}\lim_{R\rightarrow 0}\mbox{Im}\left(\int_{-R}^{-\epsilon}\frac{e^{iz}}{z}dz+\int_\epsilon^R\frac{e^{iz}}{z}dz\right)=\frac\pi2\ .[[/math]]

Now notice that

[[math]] \begin{align*} \frac1\pi\int_0^\infty \frac{e^{ius}\phi_T(-u)-e^{-ius}\phi_T(u)}{iu}du&=\frac1\pi\int_0^\infty\int_{-\infty}^\infty \frac{e^{-iu(v-s)}-e^{iu(v-s)}}{iu}dF(v)du\\ &=-\frac2\pi\int_0^\infty\int_{-\infty}^\infty \frac{\sin(u(v-s))}{u}dF(v)du\\ &=-\int_{-\infty}^\infty d F(v)\frac2\pi\int_0^\infty \frac{\sin(u(v-s))}{u}du\\ &=\int_{-\infty}^s dF(v)-\int_s^\infty dF(v)\\ &=2F(s)-1\ , \end{align*} [[/math]]

which proves \eqref{eq:gilPelaezInversion}.

■

Notice that \eqref{eq:gilPelaezInversion} can be further simplified to have

[[math]]F(s) = \frac12-\frac{1}{\pi}\int_0^\infty \Real\left[\frac{e^{-ius}\phi_T(u)}{iu}\right]du\ ,[[/math]]

which is used in pricing a European call option, as the risk-neutral probability of finishing in the money is

[[math]]\Pi_2=\mathbb Q(S_T\geq K)=\frac 12+\frac{1}{\pi}\int_0^\infty\Real \left[\frac{e^{-iu\log(K)}\phi_T(u)}{iu}\right]du[[/math]]

where [math]\mathbb Q[/math] denotes probability under the EMM. The ‘delta’ of the option is

[[math]]\Pi_1=\frac{e^{-rT}}{S_0}\mathbb E^QS_T\indicator{S_T\geq K}=\frac{e^{-rT}}{S_0}\int_{\log K}^\infty e^sf(s)ds\ ,[[/math]]

but we can define a change of measure,

[[math]]\frac{d\widetilde{\mathbb Q}}{d\mathbb Q} = \frac{e^{s_T}}{\mathbb E^Qe^{s_T}}=\frac{e^{s_T}}{S_0e^{rT}}\ ,[[/math]]

where [math]s_T=\log S_T[/math]. Then,

[[math]]\mathbb E^{\widetilde Q}e^{ius_T} = \frac{\mathbb E^ Qe^{(1+iu)s_T} }{\mathbb E^Qe^{s_T}}= \frac{\phi_T(u-i)}{\phi_T(-i)}\ ,[[/math]]

and therefore we have the formula

[[math]]\Pi_1=\widetilde{\mathbb Q}(S_T\geq K)=\frac 12+\frac{1}{\pi}\int_0^\infty \Real \left[\frac{e^{-iu\log(K)}\phi_T(u-i)}{iu\phi_T(-i)}\right]du\ .[[/math]]

Assuming no dividends, the value of a European call option with strike [math]K[/math], at time [math]t=0[/math], is

[[math]] \begin{equation} \label{eq:fourierPrice}C= S_0\Pi_1-Ke^{-rT}\Pi_2\ . \end{equation} [[/math]]

The pricing formula of equation \eqref{eq:fourierPrice} is (in one way or another) interpreted as an inverse Fourier transform. Indeed, the method of Carr and Madan shows that \eqref{eq:fourierPrice} can be well-approximated as a fast Fourier transform, which has the advantage of taking less time to compute than numerical integration, but it is also less accurate. Equation \eqref{eq:fourierPrice} has become an important piece of machinery in quantitative finance because it is (or is considered to be) an explicit solution to the pricing equation. In general the pricing formula of \eqref{eq:fourierPrice} applies to just about any situation where the characteristic function is given.

Heston Explicit Formula

Finally, an important application of the Fourier pricing formula in \eqref{eq:fourierPrice} is to the Heston model. The call option formula for the Heston model from equations \eqref{eq:dSheston} and \eqref{eq:dXheston} has a well-known inverse Fourier transform (see ) with

[[math]] \begin{align*} \Pi_j^{heston}&=\frac 12+\frac{1}{\pi}\int_0^\infty\Real \left[\frac{e^{-iu\log(K)}\phi_T^j(u)}{iu}\right]du\qquad\hbox{for }j=1,2,\\ \phi_T^j(u)&=\exp\left(A_j(u)+B_j(u)X_0+iu(\log(S_0)+rT)\right)\\ A_j(u)&=\frac{\kappa\bar X}{\gamma^2}\left((b_j-\rho \gamma ui-d_j)T-2\log\left(\frac{1-g_je^{-d_jr}}{1-g_j}\right)\right)\\ B_j(u) &= \frac{b_j-\rho \gamma ui-d_j}{\gamma^2}\left(\frac{1-e^{-d_jr}}{1-g_je^{-d_jr}}\right)\\ g_j&= \frac{b_j-\rho \gamma ui-d_j}{b_j-\rho \gamma ui+d_j}\\ d_j&= \sqrt{(\rho\gamma u i-b_j)^2-\gamma^2(2c_jui-u^2)}\\ c_1&=\frac 12,~~c_2=-\frac 12,~~b_1 = \kappa-\rho\gamma,~~b_2 = \kappa\ , \end{align*} [[/math]]

which is applied for the Heston model the volatility price of risk absorbed into parameters [math]\kappa[/math] and [math]\bar X[/math]. Warning: Use a well-made Gaussian quadrature function for numerical computation of the integrals for the Heston price; do not attempt to integrate with an evenly-spaced grid.

The Merton Jump Diffusion

Short time-to-maturity options are known to be well fit by the class of models where log-asset prices are L\'evy processes. For instance, the Merton jump diffusion,

[[math]]\frac{dS_t}{S_t} = \left(r-\nu\right)dt+\sigma dW_t+\left(e^{J_t}-1\right)dN_t[[/math]]

where [math]\sigma \gt 0[/math], [math]N_t[/math] is an independent Poisson process with parameter [math]\lambda\in(0,\infty)[/math], [math]J_t[/math] is a sequence of i.i.d. normal random variables with [math]J_{t_i}\sim N(\mu_J,\sigma_J^2)[/math] when [math]N_{t_i}-N_{t_i-} = 1[/math], and [math]\nu = \lambda\left(e^{\mu_J+\frac{\sigma_J^2}{2}}-1\right)[/math] is a compensator to ensure the process is a discounted martingale. The log-price is

[[math]] \begin{align*} \log(S_t/S_0) &= \int_0^t\left(r-\nu-\frac{\sigma^2}{2}\right)du+\sigma W_t+\int_0^tJ_udN_u\\ &=\int_0^t\left(r-\nu-\frac{\sigma^2}{2}\right)du+\sigma W_t+\sum_{n=1}^{N_t}J_n\ . \end{align*} [[/math]]

For [math]S_0=1[/math], the characteristic function is derived as follows,

[[math]] \begin{align*} \phi_T(u)&= \mathbb E\exp\left( iu\left(r-\nu-\frac{\sigma^2}{2}\right)T+iu\sigma W_T+iu\int_0^tJ_tdN_t\right)\\ &=\exp\left( iu\left(r-\nu-\frac{\sigma^2}{2}\right)T-\frac{\sigma^2u^2}{2}T+\lambda T\left(e^{iu\mu_J-\frac{\sigma_J^2u^2}{2}}-1\right)\right)\ . \end{align*} [[/math]]

This is an example of a characteristic function that is known explicitly. In addition, it is particular example of the general L\'evy-Khintchine representation,

[[math]] \begin{align*} &\mathbb E\exp\left(iu\log(S_T)\right)\\ & = \exp\left( iu\left(r-\nu-\frac{\sigma^2}{2}\right)T-\frac{\sigma^2u^2}{2}T +T\int_{\mathbb R\setminus\{0\}}\left(e^{iux}-1-iux\mathbf 1_{|x| \lt 1}\right)\eta(dx)\right)\ , \end{align*} [[/math]]

where [math]\eta[/math] is the intensity measure, satisfying [math]\int_{\mathbb R\setminus\{0\}}1\wedge x^2\eta(dx) \lt \infty[/math]. For more on the L\'evy-Khintchine representation, see .

General references

Papanicolaou, Andrew (2015). "Introduction to Stochastic Differential Equations (SDEs) for Finance". arXiv:1504.05309 [q-fin.MF].