Volatility Surface

Vol Surface: Term Structure

Term Structure usually refers to ATM implied vols by Time to maturity. And of course, these vols are all annualised vols for consistency.

Term Structure is usually upward sloping. But front-end vols are more sensitive to changes of realised vols and anticipated events (e.g. French election etc.). So market turmoil could lead to Term Structure inversion.

The Term Structure shape tend to be mean-reverting in nature. Trading strategies exploiting this mean-reverting feature involves buying and selling two options in Vega-neutral amounts, so that we have exposure to only the vols curve shape, not the level shift. This is similar to Duration-neutral Yield curve trades in Fixed Income.

Front-end vols are primarily Gamma plays, so views on Gamma is essential to formulating Term Structure. ?

Back-end vols are usually considered as sum of Front-end vols and vols curve. ?

Vol Surface: Vol Skew (Risk Reversal)

By FX market convention, Risk Reversal is quoted as

    \[\sigma _{Call, 10d} - \sigma _{Put, 10d}\]

    \[\sigma _{Call, 25d} - \sigma _{Put, 25d}\]

Risk Reversal represents directional variation of implied vol with Strike. This corresponds to the Third Standardised Central Moment of underlying spot distribution.

For reference, the n-th moment of Probability Density Function f(x) about value c is defined as:

    \[\mu _{n} = \int _{-\infty}^{\infty} (x-c)^{n}f(x) dx\]

Risk Reversal is well correlated to the correlation between spot and vol moves. We can think of Risk Reversal as the implied skew while spot-vol correlation as realised skew. And this means:

Positive skew: Option market expects spot rallies to be more volatile than sell-offs. E.g. USD/EM pairs. The option-implied spot distribution is tilted to the right.

Negative skew: Option market expects spot sell-offs to be more volatile than rallies. E.g. JPY cross pairs. The option-implied spot distribution is tilted to the left.

Skews are also often to be valued by comparing to ATM vols, i.e. \frac{RR}{ATM vol} ratio.

Vol Surface: Vol Fly (Butterfly)

By FX market convention, Butterfly is quoted as

    \[\frac{1}{2}(\sigma _{Call, 10d} + \sigma _{Put, 10d}) - \sigma _{ATM}\]

    \[\frac{1}{2}(\sigma _{Call, 25d} + \sigma _{Put, 25d}) - \sigma _{ATM}\]

Butterfly represents undirectional variation of implied vol with Strike or convexity of vol curve / smile. This corresponds to the Fourth Standardised Central Moment of underlying spot distribution.

We can think of Butterfly as the dimension of vol curve / smile that richens “wings” or low-delta options compared to ATM options. So non-Zero Butterfly means underlying spot distribution deviates from log normality assumed by Black-Scholes, and wing / low-delta vols are priced at premium to ATM vol.

Butterfly is well correlated with volatility of ATM vol. Thus, it is often considered as the parameter capturing vol-of-vol.

Vol Surface Arbitrage


Vol Surface Interpolation

Option Greeks

European Vanilla Option Greeks



  1. Option price sensitivity to spot (Spot Delta)

gradient of option price tangent line

  1. Proxy for Probability of option finishing ITM.

Delta value is between 0 and 1. But this is just a proxy for exercising probability for interpretation purpose.

In fact, the risk neutral Probability Density Function (PDF) is the 2nd derivative of call price with respect to Strike, i.e.

    \[\frac{\partial ^{2}C}{\partial K^{2}} = e^{-r(T-t)}\pi (K)\]

  1. Hedge ratio

Delta used on Vol Smile

Conventionally in FX Option space, the x-axis of Vol Smile plot is denoted as Delta (10d, 25d) instead of strikes.

Adapted Delta

Adapted Delta is the “real delta”, i.e. the actual hedge ratio, taking into account the shape of the vol smile.

Black-Scholes assumes constant vol.

Spot Delta vs Forward Delta

We normally refers Delta as Spot Delta, i.e. spot sensitivity \frac{\partial C}{\partial S}.

Forward Delta is the sensitivity to Forward price, i.e. \frac{\partial C}{\partial F} which captures interest rate risk implicit in forward points.

Forward Delta is typically used for NDF currencies and long-dated options.

Impact of Spot

Impact of Time to expiry

As we see from the below plot, as time passes, the option price curve moves closer to the At-Expiry payoff. Therefore, ITM Delta moves closer to 1 and OTM Delta moves closer to 0. ATM Delta has greater uncertainty (high Gamma) near expiry.

Impact of Vol

For ITM options, higher vol means less certainty that it will finish ITM, i.e. smaller Delta.

For OTM options, higher vol means higher probability that it will finish ITM, i.e. higher Delta.

A doubling of vol has roughly the same effect on an option’s Delta (and its price) as a quadrupling of time. For example,

    \[Call\, Option\, with\, S = 102, K = 100, T = 1m, \sigma = 5\%, Delta = 0.92\]

    \[ \sigma = 5\% \rightarrow 10\% \Rightarrow Delta = 0.92 \rightarrow 0.76\]

    \[ T = 1m \rightarrow 4m \Rightarrow Delta = 0.92 \rightarrow 0.76\]


Impact of Time to expiry

Impact of Vol

A doubling of vol has roughly the same effect on an option’s Gamma as a quadrupling of time. For example,

    \[Call\, Option\, with\, S = 100, K = 100, T = 1w, \sigma = 5\%, Delta = 0.42\]

    \[ \sigma = 5\% \rightarrow 10\% \Rightarrow Delta = 0.42 \rightarrow 0.27\]

    \[ T = 1w \rightarrow 1m \Rightarrow Delta = 0.42 \rightarrow 0.27\]

Gamma Trading

If we long an Option with Delta hedged, we will have positive P/L from long Gamma.

However, this Gamma P/L comes at cost of Theta decay as we are long option.


Theta measures the Option value decay as time passes. We say Theta is positive, meaning, as time passes (time to expiry decreases), Option price also decreases.

But can European Put Option Theta be negative?

Impact of Time to expiry

Impact of Vol


Impact of Time to expiry

Impact of Vol

Continuous Barrier Option Greeks

Reversed Knock-Out (RKO): Up and Out Call with Barrier > Strike OR Down and Out Put with Barrier < Strike

Reversed Knock-In (RKI): Up and In Call with Barrier > Strike OR Down and In Put with Barrier < Strike

Reversed means the Barrier level is In-The-Money (ITM).

RKO Call / Put + RKI Call / Put = Vanilla Call / Put

This is true for option price and all greeks.


Impact of Spot

Delta is the gradient of curve. We notice there’s a Delta gap on barrier trigger, for both RKI and RKO.

RKI Call Option

RKO Call Option

American Vanilla Option Pricing – Binomial Tree Method

We first divide the American Call option tenor into smaller time steps, each represented as \Delta t.

In each time step, assume underlying asset price may move from initial value S either up to Su with real-world probability q or down to Sd with 1-q.

We assume the annualised risk free rate is r.

At the end of this time step \Delta t, the payoff of the American Call option is

    \[C_{u} = Max(0, Su-K), if\, underlying\, goes\, up\, to\, Su\]

    \[C_{d} = Max(0, Sd-K), if\, underlying\, goes\, down\, to\, Sd\]

We then construct a portfolio with h shares of underlying asset and D amount of cash invested at risk free rate. The initial portfolio cost is hS + D, the portfolio value at end of the time step is hSu + De^{r\Delta t} or hSd + De^{r\Delta t} depending on underlying.

And we can carefully choose below h and D so that this portfolio replicates the payoff of American Call option at end of this time step.

    \[h = \frac{C_{u} - C_{d}}{(c - d)S}\]

    \[D = \frac{uC_{d} -dC_{u}}{(u - d)e^{r\Delta t}}\]

With no arbitrage, the value of the American Call option has to be equal to the portfolio initial cost at beginning of the time step. So we have:

    \[C = hS + D = \frac{C_{u} - C_{d}}{(c - d)} + \frac{uC_{d} -dC_{u}}{(u - d)e^{r\Delta t}} = \frac{pC_{u} + (1-p)C_{d}}{e^{r\Delta t}}\]


    \[p = \frac{e^{r\Delta t} - d}{u - d}\]

    \[1 - p = \frac{u - e^{r\Delta t}}{u - d}\]

So we notice that the real-world probability q is NOT in the formula, which means the American Call option prices does not depend on investors’ individual risk preference.

We also notice p is between 0 and 1, and therefore be regarded as the risk neutral probability. And the American Call option price is just the discounted value of future expected payoffs.

In the risk neutral world, p is the probability that underlying asset price goes up to Su, and 1-p is the probability it goes down to Sd. And the current underlying price is just the discounted value of future expected payoffs, i.e.

    \[S = \frac{pSu + (1-p)Sd}{e^{r\Delta t}}\]

which can be simplified as:

    \[e^{r\Delta t} = pu + (1-p)d\]

From our assumption, we know the underlying asset price return over one time step \Delta t follows Binomial distribution, thus the one-step variance of price return is \sigma ^{2}\Delta t.

    \[pu^{2}+(1-p)d^{2} - [pu+(1-p)d]^{2} = \sigma ^{2}\Delta t\]

Cox introduced another condition:

    \[u = \frac{1}{d}\]

Now, from the above 3 equations, we can solve for p,\, u,\, d for given r,\, \sigma,\, \Delta t as below:

    \[u = e^{\sigma \sqrt{t}}\]

    \[d = e^{-\sigma \sqrt{t}}\]

    \[p = \frac{e^{r\Delta t} - d}{u - d}\]

So now we can use above Binomial Tree model to calculate the American Call option price step by step backward from expiry. At each step, we need to evaluate the discounted future expected payoffs and the payoff if we exercise at this step. We take the bigger value of the two as the value for this time step node.

Python code:

def binomialTree(callPut, spot, strike, rate, sigma, tenor, N=2000, american=True):
    # Each time step period
    deltaT = float(tenor) / N
    u = np.exp(sigma * np.sqrt(deltaT))
    d = 1.0 / u
    a = np.exp(rate * deltaT)
    p = (a - d) / (u - d)
    oneMinusP = 1.0 - p
    # Initialize the arrays
    fs = np.asarray([0.0 for i in xrange(N + 1)])
    # Stock tree for calculations of expiration values
    fs2 = np.asarray([(spot * u ** j * d ** (N - j)) for j in xrange(N + 1)])
    # Vectorize the strikes to speed up expiration check
    fs3 = np.asarray([float(strike) for i in xrange(N + 1)])
    # Compute the Binomial Tree leaves, f_{N, j}
    if callPut == 'Call':
        fs[:] = np.maximum(fs2 - fs3, 0.0)
        fs[:] = np.maximum(-fs2 + fs3, 0.0)
    # Calculate backward the option prices
    for i in xrange(N - 1, -1, -1):
        fs[:-1] = np.exp(-rate * deltaT) * (p * fs[1:] + oneMinusP * fs[:-1])
        fs2[:] = fs2[:] * u
        if american:
            # Simply check if the option is worth more alive or dead
            if callPut == 'Call':
                fs[:] = np.maximum(fs[:], fs2[:] - fs3[:])
                fs[:] = np.maximum(fs[:], -fs2[:] + fs3[:])
    return fs[0]

European Vanilla Option Pricing – Black-Scholes PDE

Assume underlying spot follows Geometric Brownian Motion, i.e.

    \[dS = \mu Sdt + \sigma SdW_{t}\]

Let C be the call option price. We obtain dC using Ito Lemma

    \[dC = [\frac{\partial C}{\partial t} + \frac{\partial C}{\partial S}\mu S + \frac{1}{2}\frac{\partial ^{2} C}{\partial S^{2}}\sigma ^{2}S^{2}]dt + \frac{\partial C}{\partial S}\sigma Sdz\]

Construct a delta neutral portfolio (short call option and long underlying), then we have:

    \[dV = \frac{\partial C}{\partial S}dS - dC\]

If we combine the terms, we will get

    \[dV = [-\frac{\partial C}{\partial t} - \frac{1}{2}\frac{\partial ^{2}C}{\partial S^{2}} \sigma ^{2}S^{2}]dt\]

Realise dV is independent of random term dz, thus portfolio V is risk free.
Realise dV is independent of expected return \mu.

Thus, portfolio V should earn the risk free rate of return, i.e.

    \[dV = rVdt = r(\frac{\partial C}{\partial S}S - C)dt\]

Therefore, combining with dV in the previous step, we have below Black-Scholes PDE:

    \[\frac{\partial C}{\partial t} + \frac{1}{2}\frac{\partial ^{2}C}{\partial S^{2}}\sigma ^{2}S^{2} + \frac{\partial C}{\partial S}rS = rC\]

Now we need to solve the above Black-Scholes PDE.

Step 1

Transformation: Let’s introduce new variables x = \ln{\frac{S}{K}}, and \tau = T-t.

Therefore, the Call option price C(S, t) can be represented using new variables x and \tau as C(Ke^{x}, T-\tau).

Now we introduce a new function Z(x, \tau) = C(Ke^{x}, T-\tau). We need to find the PDE for Z(x, \tau) where x \in \mathbb{R}, \tau \in [0, T]

By Chain rule for partial derivatives, we have:

    \[\frac{\partial C}{\partial S} = \frac{\partial Z}{\partial x} \frac{\partial x}{\partial S} + \frac{\partial Z}{\partial \tau} \frac{\partial \tau}{\partial S}=\frac{\partial Z}{\partial x} \frac{K}{S} \frac{1}{K}=\frac{\partial Z}{\partial x} \frac{1}{S}\]

    \[\frac{\partial ^{2} C}{\partial S^{2}} = \frac{\partial (\frac{\partial Z}{\partial x} \frac{1}{S})}{\partial S} = \frac{\partial (\frac{\partial Z}{\partial x}) \frac{1}{S}}{\partial S} + \frac{\frac{\partial Z}{\partial x} \partial (\frac{1}{S})}{\partial S} \]

    \[= [\frac{\partial (\frac{\partial Z}{\partial x})}{\partial x} \frac{\partial x}{\partial S} + \frac{\partial (\frac{\partial Z}{\partial x})}{\partial \tau} \frac{\partial \tau}{\partial S}] \frac{1}{S} + \frac{\partial Z}{\partial x}(-\frac{1}{S^{2}}) = \frac{\partial ^{2}Z}{\partial x^{2}} \frac{1}{S^{2}} - \frac{\partial Z}{\partial x} \frac{1}{S^{2}}\]

    \[\frac{\partial C}{\partial t} = \frac{\partial C}{\partial \tau}\frac{\partial \tau}{\partial t} + \frac{\partial Z}{\partial x} \frac{\partial x}{\partial t} = - \frac{\partial Z}{\partial \tau}\]

Now we plug \frac{\partial C}{\partial S}, \frac{\partial ^{2} C}{\partial S^{2}}, \frac{\partial C}{\partial t}, Z into the Black-Scholes PDE, then we find the PDE for Z(x, \tau):

    \[\frac{\sigma ^{2}}{2} \frac{\partial ^{2} Z}{\partial x^{2}} + (r-\frac{\sigma ^{2}}{2})\frac{\partial Z}{\partial x} -\frac{\partial Z}{\partial \tau} - Zr = 0\]

Step 2

Transformation to Heat Equation: Let’s introduce a new function u(x, \tau) = e^{\alpha x + \beta \tau}Z(x,\tau). We need to choose constants \alpha, \beta \in \mathbb{R} so that the PDE of u is Heat Equation.

    \[\frac{\partial u}{\partial \tau} = \beta e^{\alpha x + \beta \tau} Z + e^{\alpha x + \beta \tau}\frac{\partial Z}{\partial \tau} \]

    \[\frac{\partial u}{\partial x} = \alpha e^{\alpha x + \beta \tau} Z + e^{\alpha x + \beta \tau}\frac{\partial Z}{\partial x} \]

    \[\frac{\partial ^{2}u}{\partial x^{2}} = \alpha ^{2} e^{\alpha x + \beta \tau}Z +\alpha e^{\alpha x + \beta \tau} \frac{\partial Z}{\partial x} + \alpha e^{\alpha x + \beta \tau} \frac{\partial Z}{\partial x} + e^{\alpha x + \beta \tau} \frac{\partial ^{2}Z}{\partial x^{2}} \]

    \[= e^{\alpha x + \beta \tau}\frac{\partial ^{2} Z}{\partial x^{2}} + 2\alpha e^{\alpha x + \beta \tau}\frac{\partial Z}{\partial x} + \alpha ^{2}e^{\alpha x + \beta \tau}Z \]

Together with the PDE for Z, we can derive the PDE for u:

    \[\frac{\partial u}{\partial \tau} -\frac{\sigma ^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}} + (\alpha \sigma ^{2} + \frac{\sigma ^{2}}{2} - r)\frac{\partial u}{\partial x} + (r + r\alpha -\frac{\sigma ^{2} \alpha ^{2}}{2} - \frac{\alpha \sigma ^{2}}{2})u = 0\]

To be a Heat Equation, we need to force the last two terms be 0. Thus

    \[\alpha \sigma ^{2} + \frac{\sigma ^{2}}{2} - r = 0\]

    \[r + r\alpha -\frac{\sigma ^{2} \alpha ^{2}}{2} - \frac{\alpha \sigma ^{2}}{2} = 0\]

Then we have

    \[\alpha = \frac{r}{\sigma ^{2}} - \frac{1}{2}\]

    \[\beta = \frac{r}{2} + \frac{\sigma ^{2}}{8} + \frac{r^{2}}{2\sigma ^{2}}\]

Step 3

The solution u(x, \tau) of PDE \frac{\partial u}{\partial \tau} -\frac{\sigma ^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}} = 0 is given by Green formula as below:

    \[u(x, \tau) = \frac{1}{\sqrt{2\sigma ^{2}\pi \tau}} \int _{-\infty}^{\infty} e^{-\frac{(x-s)^{2}}{2\sigma ^{2} \tau}}u(s, 0) ds \]

Step 4

We look at the boundary condition u(x, 0).

    \[u(x, 0) = e^{\alpha x}Z(x, 0) = e^{\alpha x}C(S, T) =  e^{\alpha x}(S-K) \,if \,x>0\]

    \[u(x, 0) = 0 \,if \,x \leq 0\]


    \[u(x, \tau) = \frac{1}{\sqrt{2\sigma ^{2}\pi \tau}} \int _{0}^{\infty} e^{-\frac{(x-s)^{2}}{2\sigma ^{2} \tau}}e^{\alpha s}(S-K) ds \]

which can be integrated as below, where \phi is the cumulative distribution function (CDF) for Normal distribution.

    \[u(x, \tau)=Se^{\alpha x + \beta \tau} \phi (\frac{x+(r+\frac{\sigma ^{2}}{2})\tau}{\sigma \sqrt{\tau}}) - Ke^{\alpha x + \frac{1}{2}\sigma ^{2}\tau \alpha ^{2}} \phi (\frac{x+\sigma ^{2}\tau \alpha}{\sigma \sqrt{\tau}}) \]

Step 5

From the above steps, we have relation

    \[u(x, \tau) = e^{\alpha x + \beta \tau}Z(x,\tau) = e^{\alpha x + \beta \tau}C(Ke^{x}, T-\tau) = e^{\alpha x + \beta \tau}C(S, t)\]

And from Step 4, we know the result of u(x, \tau).

Therefore, we derive C(S, t) as

    \[C(S, t) = S\phi (\frac{x+(r+\frac{\sigma ^{2}}{2})\tau}{\sigma \sqrt{\tau}}) - Ke^{\frac{1}{2}\sigma ^{2}\tau \alpha ^{2}-\beta \tau} \phi (\frac{x+\sigma ^{2}\tau \alpha}{\sigma \sqrt{\tau}}) \]

Now we plug in x, \tau, \alpha, \beta from previous steps. Finally, Call option price C(S, t) can be represented as

    \[C(S, t) = S\phi (d_{1}) - e^{-r(T-t)}K\phi (d_{2})\]


    \[d_{1} = \frac{\ln \frac{S}{K} + (r+\frac{\sigma ^{2}}{2})(T-t)}{\sigma \sqrt{T-t}} \]

    \[d_{2} = d_{1} - \sigma \sqrt{T-t} = \frac{\ln \frac{S}{K} + (r-\frac{\sigma ^{2}}{2})(T-t)}{\sigma \sqrt{T-t}}\]


Python implementation of Black-Scholes formula:

def ncdf(x):
    Cumulative distribution function for the standard normal distribution.
    Alternatively, we can use below:
    from scipy.stats import norm
    return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
def npdf(x):
    Probability distribution function for the standard normal distribution.
    Alternatively, we can use below:
    from scipy.stats import norm
    return np.exp(-np.square(x) / 2) / np.sqrt(2 * np.pi)
def blackScholesOptionPrice(callPut, spot, strike, tenor, rate, sigma):
    Black-Scholes option pricing
    tenor is float in years. e.g. tenor for 6 month is 0.5
    d1 = (np.log(spot / strike) + (rate + 0.5 * sigma ** 2) * tenor) / (sigma * np.sqrt(tenor))
    d2 = d1 - sigma * np.sqrt(tenor)
    if callPut == 'Call':
        return spot * ncdf(d1) - strike * np.exp(-rate * tenor) * ncdf(d2)
    elif callPut == 'Put':
        return -spot * ncdf(-d1) + strike * np.exp(-rate * tenor) * ncdf(-d2)
def blackScholesVega(callPut, spot, strike, tenor, rate, sigma):
    """ Black-Scholes vega """
    d1 = (np.log(spot / strike) + (rate + 0.5 * sigma ** 2) * tenor) / (sigma * np.sqrt(tenor))
    return spot * np.sqrt(tenor) * npdf(d1)
def blackScholesDelta(callPut, spot, strike, tenor, rate, sigma):
    """ Black-Scholes delta """
    d1 = (np.log(spot / strike) + (rate + 0.5 * sigma ** 2) * tenor) / (sigma * np.sqrt(tenor))
    if callPut == 'Call':
        return ncdf(d1)
    elif callPut == 'Put':
        return ncdf(d1) - 1
def blackScholesGamma(callPut, spot, strike, tenor, rate, sigma):
    """" Black-Scholes gamma """
    d1 = (np.log(spot / strike) + (rate + 0.5 * sigma ** 2) * tenor) / (sigma * np.sqrt(tenor))
    return npdf(d1) / (spot * sigma * np.sqrt(tenor))


进入2018年2月以来市场打破了几年来的平静,有据说程序带来的闪崩,有FB为主的data privacy的担心,更有特朗普炮轰亚马逊和trade war的恐慌。总之,一个字:不平静。


比如,如果我现在已经massive long FB,但它的股票一直在跌,我现在就不可以long call 或者short put,而是应该相反。hedge的size和tenor也要根据自己的view来决定。


尤其在time horizon很长(比如10年)和leverage很低的情况下。这种情况下不需要担心爆仓等。如果是passive investment,比如buy and hold S&P500,一定要相信经济长期是盘旋向上的,尤其美国经济也是整体向上的,这在短时间内不会改变。即便是明年经济危机,最多两年又会恢复。我们从Cobb-Douglas production function就可以看出productivity, capital, labor都会推动整个经济增长,而股市往往又是跟随经济起伏。记得之前有人十年如一日做空VIX发家致富大概也是一个道理。

Perfect hedging

这是另外一种极端的情况。比如对于卖方/做市商来讲,他们通常只赚取spread和flow,不会take risk。假如一个FX Derivatives market maker有residual risk exposure (delta, gamma, vega等等),他会马上去做hedging,让自己的PNL不会因spot, vol波动而起伏。


  1. 假如我现在long 10K USD SPY spot position,我的view是2018年会比较波动但不会有大危机,预期5%-10%的回报率。这种情况下,我就可以做一个collar,short out-of-money call 同时 long out-of-money put,这样的hedge不会很贵因为premium抵消一些,可以做到获取一定范围内的upsides收益,同时也有downside protection.
  2. 假如我的view是美股今年可能继续大涨,downside risk不会很大。那我可以在现有Long spot position基础上,Long higher strike OTM put, 同时short lower strike OTM put. 这样我们去除了一部分downside risk,同时hedge不会很贵。唯一的风险是美股大幅下跌的情况,但根据我们的view是小概率事件。
  3. 假如我现在Long FB股票,但市场对科技股情绪下降,如果我相信FB的前景和基本面好于其他科技股,如Twitter,特斯拉等,我可以选择short 这些科技股来对冲掉科技股行业的风险,只留下company specific risk。这就是Long short的思想。或者我可以short tech index,因为我相信FB会outperform 行业。



Short option一定要非常谨慎

因为有unlimited downside risk。而我们往往低估或不能预测市场的走向和不确定性。08年经济危机的时候就有人利用这点赚了很多钱。危机之所以是危机,是因为它出乎人们的预料。

European Vanilla Option Pricing – Monte Carlo Methods

  1. We start with the assumption that underlying follow Geometric Brownian Motion (GBM):

    \[dS_{t} = \mu S_{t}dt + \sigma S_{t}dW_{t}\]

  1. We use Ito’s Lemma with G = \ln (S), then we have

    \[\frac{\partial G}{\partial S} = \frac{1}{S}\]

    \[\frac{\partial ^{2}G}{\partial S^{2}} = -\frac{1}{S^{2}}\]

    \[\frac{\partial G}{\partial t} = 0\]

By Ito’s Lemma, we have

    \[dG = (\mu - \frac{\sigma ^{2}}{2})dt + \sigma dz\]

  1. Therefore, the change of \ln(S) between time 0 and future time T, is normally distributed as following:

    \[\ln(S_{T}) - \ln(S_{0}) \sim \mathcal{N}[(\mu - \frac{\sigma ^{2}}{2})T, \sigma \sqrt{T}]\]

Thus, the future underlying price can be written as,

    \[S_{T} = S_{0}e^{(\mu-\frac{\sigma ^{2}}{2})T + \sigma \sqrt{T}\epsilon}\]

\epsilon is the noise term from standard normal distribution.

Note, we will take \mu = r, which is the risk free rate. This means investors are risk neutral and requires risk free return on underlying asset. This is to be consistent with the risk neutral probabilities used in simulation.
Correspondingly, we also use risk free rate in the discount factor in step 5.

  1. So now we can simulate the future underlying price at expiry. With European Call or Put boundary condition to calculate the payoff.

    \[Call Payoff at Expiry = Max[(S_{T} - K), 0]\]

    \[Put Payoff at Expiry = Max[(K - S_{T}), 0]\]

  1. We then need to discount the future payoff back to present by multiplying a discount factor,

    \[Discount Factor = e^{-rT}\]

  1. The above two steps are repeated many times and its expectation is calculated as the final simulation result.
python code:
def getMCPrice(self):
    Determine the option price using a Monte Carlo approach.
    The log return of underlying follow Normal distribution.
    s_T = s_t * exp((r - 1/2 * sig^2) * (T-t) + sig * sqrt(T-t) * sig_Normal)
    calc = np.zeros([self.iterations, 2])
    rand = np.random.normal(0, 1, [1, self.iterations])
    mult = self.spot * np.exp(self.tenor * (self.rate - 0.5 * self.sigma**2))
    if self.callPut == 'Call':
        calc[:,1] = mult * np.exp(np.sqrt((self.sigma**2)*self.tenor) * rand) - self.strike
    elif self.callPut == 'Put':
        calc[:,1] = self.strike - mult*np.exp(np.sqrt((self.sigma**2) * self.tenor) * rand)
    avgPayOff = np.sum(np.amax(calc, axis=1)) / float(self.iterations)
    return np.exp(-self.rate * self.tenor) * avgPayOff
def getBSPrice(self):
    """ Determine the option price using the exact Black-Scholes expression. """
    return blackScholesOptionPrice(self.callPut, self.spot, self.strike, self.tenor, self.rate, self.sigma)
We can run the above in Python console:
from option_pricer import EuropeanVanillaPricer
pricer = EuropeanVanillaPricer()
As we can see in the above Monte Carlo simulation, we rely on drawing random numbers, \epsilon from a Standard Normal distribution.
Alternatively, we can use random numbers from a Uniform distribution, i.e. equal probability of each random number.


To do this, we combine step 3, 4 and 5, the current option price is obtained by integrating the terminal payoff under the risk neutral measure:

    \[e^{-rT}\int_{-\infty}^{\infty}g(S_{0}e^{(\mu-\frac{\sigma ^{2}}{2})T + \sigma \sqrt{T} \epsilon}) \frac{1}{\sqrt{2\pi}}e^{-\frac{\epsilon ^{2}}{2}}d\epsilon\]

    \[=\int_{-\infty}^{\infty}h(\epsilon)\frac{1}{\sqrt{2\pi}}e^{-\frac{\epsilon ^{2}}{2}}d\epsilon\]

    \[=\int_{0}^{1}h[\phi ^{-1}(x)]dx\]


In the first line, function g is just the payoff condition at expiry. As we are integrating with regard to \epsilon, which follows Standard Normal distribution, the last term is the probability density function.
In the second line, we just use a new function h of epsilon to make the expression more compact.
In the third line, we do inverse transformation to integrate with regard to the cumulative probability, x.


So now it becomes an integral of function f over the UNIFORM distribution with range [0, 1].
Now our simulation task becomes taking random number x from the Uniform distribution [0, 1], and then calculate integral of function f using Monte Carlo.


To be more specific, our task has been changed from calculating

    \[\frac{1}{N}\sum_{i=1}^{N}h(\epsilon _{i}), where \, \epsilon _{i} \sim \mathcal{N}(0, 1)\]

To evaluating



When we use Monte Carlo to estimate function integral, we may run into problem of random number clustering, which essentially leads to Convergence rate of O(N^{0.5}).
To conquer this issue, instead of using pseudo-random numbers, we can use a deterministic sequence, whose numbers are more equally spaced. And this is exactly what Quasi Monte Carlo (QMC) does. More details on the low-discrepancy sequences can be found in this post.


class VolatilityPricer():
    Realized vol:
    Same as Black-Scholes, we assume the underlying follows a Geometric Brownian Motion.
    Then its log return follows a Normal distribution, with mean as 0.
    We take as input the historical daily underlying prices.
    Annualization factor is 252.
    Degree of Freedom is 0 as we are calculating the exact realized vol for the given historical period.
    Implied vol:
    Use Black-Scholes to back out the implied volatility from the given market option price.
    def __init__(self):
        self.historicalDataBySymbol = dict()
        self.dataHub = DataHub()
        self.realizedVolBySymbol = dict()
    def _loadHistoricalUnderlyingData(self, startDate, endDate, symbols):
        self.historicalDataBySymbol = self.dataHub.downloadDataFromYahoo(startDate, endDate, symbols)
    def _calculateRealizedVol(self, ts):
        """ Calculate the realized vol from given time series """
        pctChange = ts.pct_change().dropna()
        logReturns = np.log(1+pctChange)
        vol = np.sqrt(np.sum(np.square(logReturns)) / logReturns.size)
        annualizedVol = vol * np.sqrt(252)
        return annualizedVol
    def getRealizedVol(self, startDate=datetime.date.today()-datetime.timedelta(days=30), endDate=datetime.date.today(), symbols=['SPY']):
        """ Calculate the realized volatility from historical market data """
        self._loadHistoricalUnderlyingData(startDate, endDate, symbols)
        for symbol, df in self.historicalDataBySymbol.iteritems():
            # Use daily Close to calculate realized vols
            realizedVol = self._calculateRealizedVol(df.loc[:, 'Close'])
            self.realizedVolBySymbol[symbol] = realizedVol
        return self.realizedVolBySymbol
    def getImpliedVol(self, optionPrice=17.5, callPut='Call', spot=586.08, strike=585.0, tenor=0.109589, rate=0.0002):
        """ Calculate the implied volatility from option market price """
        return blackScholesSolveImpliedVol(optionPrice, callPut, spot, strike, tenor, rate)


这个我们也是假设股票的价格遵循Geometric Browian Motion,进而它的log return就服从Normal Distribution。
我的pricer可以自动从Yahoo Finance抓取历史数据,并以此构建时间序列。这样我们就可以算出股票的log returns序列,进而就可以算出它的方差,最后把波动率年化。
from volatility_pricer import VolatilityPricer
vp = VolatilityPricer()
{'SPY': 0.086197389793546381}
{'SPY': 0.16562165494524139}


这个原理就很简单了,用Black-Scholes已知期权价格back out波动率。
from volatility_pricer import VolatilityPricer
vp = VolatilityPricer()