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}}\]

with

    \[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)
    else:
        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[:])
            else:
                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\]

Then

    \[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})\]

where

    \[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
    norm.cdf(x)
    """
    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
    norm.pdf(x)
    """
    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))

Monte Carlo Methods

Why

Simulation / Compute function expectation

Compute function integral

Suppose we want to find the value of:

    \[\int_{a}^{b}f(x)dx\]

in some region with volume V.
Now, we can estimate this integral by estimating the proportion of random points that fall under f(x), then multiplied by V.

Important in Exotics pricing as there’s no analytical solutions (e.g. Asian Options)

Convergence Rate

The speed Monte Carlo method converges to the correct result as we increase the sample size / number of simulations.

Monte Carlo has convergence rate O(N^{0.5}).

Quasi Monte Carlo has convergence rate O(\frac{1}{N}).

Variance Reduction

Reduce the variance of Monte Carlo simulated result with regard to the true result.

Increase number of simulations

Not very feasible, as in order to decrease variance / noise linearly, we have to increase simulations exponentially.

Importance Sampling

 

Quasi Monte Carlo (QMC)

Low-discrepancy Sequence
Halton Sequence
Sobol Sequence
Faure Sequence

We take any prime number r where r>=2. Any integer n has a unique expansion with base r. We can then generate a sequence of numbers in the interval [0, 1), which are equally spaced within the interval.

For example, r = 3 and n = 7. We can write 7 in the form of base 3 as below:

    \[7 = 2(3^{1}) + 1(3^{0}) = 21_{3}\]

Now, if we reflect this number about its “decimal point”, we get a new number in [0, 1):

    \[1(3^{-1}) + 2(3^{-2}) = \frac{5}{9}\]

We keep on doing this for every number n, we will generate a sequence in interval [0, 1). And we observe that the newly generated number keep filling the gaps in proceeding sequence. For example, for n = 1 to 9, we have below Faure Sequence:

    \[\frac{9}{27}, \frac{18}{27}, \frac{3}{27}, \frac{12}{27}, \frac{21}{27}, \frac{6}{27}, \frac{15}{27}, \frac{24}{27}, \frac{1}{27}, \]

In a general form, for any given number n, we can represent n in base r:

    \[n = \sum_{j=0}^{m}a_{j}(n)r^{j}\]

Then, we can find the corresponding Faure number in interval [0, 1) as follows:

    \[\sum_{j=0}^{m}a_{j}(n)r^{-j-1}\]

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()
pricer.getMCPrice()
2.1620364099067015
pricer.getBSPrice()
2.1736062697012564
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\]

    \[=\int_{0}^{1}f(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

    \[=\int_{0}^{1}f(x)dx\]

 

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.

用Python计算分析实现波动率和隐含波动率

今天又敲了一个volatility_pricer.py,可以分别计算给定股票的实现波动率,和用Black-Scholes算出的隐含波动率。
Python代码可以详见:
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()
vp.getRealizedVol()
{'SPY': 0.086197389793546381}
vp.getRealizedVol(startDate=datetime.date(2018,1,1))
{'SPY': 0.16562165494524139}

计算给定期权的隐含波动率

这个原理就很简单了,用Black-Scholes已知期权价格back out波动率。
但是没有公式可以直接进行计算,所以我用了牛顿导数方法来不断求最优解。我们在每一步需要计算一个当前波动率下的期权价格以及对应的Vega,然后不断缩小误差,直到满意的误差范围。
使用方法:
from volatility_pricer import VolatilityPricer
vp = VolatilityPricer()
vp.getImpliedVol()
0.21921387741959775
哈哈,我还是一如既往地太懒写细节,大家就自己看code吧。
如果有问题或者需要交流,欢迎关注和联系我哈。