Put-Call Parity

Put-Call parity is a simple result connecting the prices of puts and calls in a model-independent way via the forward price.

Consider the three graphs below, showing independently the payoff at expiry of a vanilla call, a vanilla put, and a forward contract. We can see that the payoff of a long call plus the payoff of a short put will precisely overlap the forward contract payoff (assuming they have the same strike and expiry…).

The combination of a long call and a short put with the same strike and expiry is equivalent to a forward at the same strike and expiry. This is guaranteed by their payoffs at expiry (making the usual assumptions about tradability of the underlying/forward etc.), since the payoff of holding a long call and a short strike can be exactly hedged by holding a forward contract
The combination of a long call and a short put with the same strike and expiry is equivalent to a forward at the same strike and expiry. This is guaranteed by their payoffs at expiry (making the usual assumptions about tradability of the underlying/forward etc.), since the payoff of holding a long call and a short strike can be exactly hedged by holding a forward contract

This can be seen from the algebra as well:

\begin{align*} C_{\rm call}(T) - C_{\rm put}(T)\ & = \big( S(T) - K \big)^+ - \big( K - S(T) \big)^+ \\ & = S(T) - K \\ & = F(T,T) \end{align*}

If two portfolios have the same payoff, it’s a fundamental rule of derivatives pricing that they must have the same price, which gives the fundamental put-call parity relationship

C_{\rm call}(t) - C_{\rm put}(t)\ = Z(t,T)\cdot\Big( F(t,T) - K\Big)

We can take the BS expressions for put and call prices given in a previous post and observe that they obey the relationship [since \inline \Phi(x) + \Phi(-x) = 1], but importantly this result is model independent. As long as there is a forward contract available with the strike and expiry of the options, then this result MUST hold in ANY model.

One implication of this is that if we use a model to determine the price of a call option, put-call parity fixes the price of the put option to be

C_{\rm put}(t) = C_{\rm call}(t) - Z(t,T)\cdot\Big( F(t,T) - K\Big)

We can actually use this to improve our code. Out-of-the-money options typically require more paths to converge in Monte Carlo simulations because they rely on the detailed behaviour of a few paths that travel a long way. However, we can calculate in-the-money prices quickly and due to put-call parity these will lock the price of the matching out-of-the-money options.

Put-call parity is also an important test of the implementation of a model – if the prices that we are getting out don’t obey this relationship then we’ve done something seriously wrong [as long as both prices have converged!]. An equivalent statement is that the implied BS volatility of matching puts and calls should be exactly the same in any model, since the implied vol is the BS vol that would produce the given price for a put/call, and this is locked by put-call parity.

We can get some other similar relationships between put and call prices that I’ll look at in the future. Also note that put-call parity holds for european-style options, but not for many more complicated path-dependent options or those with early exercise features like American options.

Stochastic Rates Models

In a previous post, I introduced the Black-Scholes model, in which the price of an underlying stock is modeled with a stochastic variable that changes unpredictably with time. I’ve also discussed the basic model-independent rates products whose value can be determined at the present time exactly. However, to progress further with interest rate derivatives, we’re going to need to model interest rates more carefully. We’ve assumed rates are deterministic so far, but of course this isn’t true – just like stocks, they change with time in an uncertain manner, so we need to allow them to become stochastic as well.

One way of doing this is by analogy with the BS case, by allowing the short rate (which is the instantaneous risk-neutral interest rate r(t)) to become stochastic as well, and then integrating it over the required period of time to calculate forward rates.

A very basic example of this is the Vasicek Model. In this model the short rate is defined to be stochastic, with behaviour governed by the following SDE

    \[dr = k\ [\ \theta - r(t)\ ]\ dt\ +\ \sigma\ dW\]

where k\theta and \sigma are constants and dW is the standard Wiener increment as described before. This is marginally more complicated than the BS model, but still belongs to the small family of SDEs that are analytically tractable. Unlike stock prices, we expect rates to be mean-reverting – stock price variance grows with time, but we expect the distribution of rates to be confined to a fairly narrow range by comparison. The term in square brackets above achieves this, since if r(t) is greater than \theta then it will be negative and cause the rate to be pulled down, while if it is below \theta the term will be positive and push the rate up towards \theta.

Solving the equation requires a trick, which is instead of thinking about the rate alone to think about the quantity d\big(\ e^{kt}\cdot r(t)\ \big). This is equal to e^{kt}\cdot dr(t)\ + \ ke^{kt} dt \cdot r(t), and substituting in the incremental rate term from the original equation we have

    \[d\big(\ e^{kt}\cdot r(t)\ \big) = e^{kt} \Big(k\ [ \theta - r(t) ] dt + \sigma dW \Big) + \ e^{kt} dt \cdot r(t)\]

    \[d\big(\ e^{kt}\cdot r(t)\ \big) = e^{kt} \Big( k\cdot \theta\cdot dt + \sigma\cdot dW \Big)\]

note that the term in r(t) has been cancelled out, and the remaining terms can be integrated directly from a starting time s to a finishing time T

    \[\Big[\ e^{kt}\cdot r(t)\ \Big]_{t=s}^{t=T} = \int_{t=s}^{t=T} k\ \theta\ e^{kt}\cdot dt + \int_{t=s}^{t=T} e^{kt} \sigma\cdot dW\]

    \[e^{kT}\cdot r(T)\ - e^{ks}\cdot r(s)\ = \Big[\ \theta\ e^{kt}\ \Big]_{t=s}^{t=T} + \sigma \int_{t=s}^{t=T} e^{kt}\cdot dW\]

    \[r(T)\ = e^{-k(T-s)}\cdot r(s)\ + \theta\Big(\ 1 - e^{-k(T-s)}\ \Big) + \sigma \int_{s}^{T} e^{-k(T-t)}\cdot dW(t)\]

where r(s) is the initial rate. This is simply a gaussian distribution with mean and variance given by

    \[\mathbb{E} \big[\ r(T)\ |\ {\cal F}_{\rm s}\ \big] = e^{-k(T-s)}\cdot r(s)\ +\ \theta\Big(\ 1 - e^{-k(T-s)}\ \Big)\]

    \[\mathbb{V} \big[\ r(T)\ |\ {\cal F}_{\rm s}\ \big] = {\sigma^2 \over 2k} \Big(\ 1 - e^{-2k(T-s)}\ \Big)\]

Where the variance is calculated using the “Ito isometry

    \[\mathbb{V} \big[\ r(T)\ |\ {\cal F}_{\rm s}\ \big] = \mathbb{E} \big[\ \big(r(T)\big)^2\ |\ {\cal F}_{\rm s}\ \big] - \big(\ \mathbb{E} \big[\ r(T)\ |\ {\cal F}_{\rm s}\ \big]\ \big)^2\]

    \[= \Big( \sigma \int_s^T e^{-k(T-t)} dW(t) \Big)^2\]

    \[=\sigma^2 \int_s^T \Big( e^{-k(T-t)}\Big)^2 dt\]

    \[= {\sigma^2 \over 2k} \Big(\ 1 - e^{-2k(T-s)}\ \Big)\]

as stated above. Note that this allows the possibility of rates going negative, which is generally considered to be a weakness of the model, but the chance is usually rather small.




As we know the distribution of the short rate, we can calculate some other relevant quantities. Of primary importance are the Zero Coupon Bonds, which are required for calculation of forward interest rates. A ZCB is a derivative that pays $1 at a future time t, and we can price this using the Risk-Neutral Valuation technique.

According to the fundamental theorem of asset pricing, the current price of a derivative divided by our choice of numeraire must be equal to it’s future expected price at any time divided by the value of the numeraire at that time, with the expectation taken in the probability measure corresponding to the choice of numeraire.

In the risk-neutral measure, the numeraire is just a unit of currency, initially worth $1 but continually re-invested at the instantaneous short rate, so that its price at time t is $1 \cdot e^{\int_0^t r(t') dt}. Now, the price of the ZCB is given by

    \[{Z(0,t) \over \$1} = {\mathbb E}^{RN}\Big[\ {Z(s,t)\over \$1 \cdot e^{\int_{0}^s r(t') dt' }}\ |\ {\cal F}_0\ \Big]\]

    \[Z(0,t) = {\mathbb E}\Big[\ {Z(s,t)\cdot e^{-\int_{0}^s r(t') dt' }}\ |\ {\cal F}_0\ \Big]\]

Although the RHS is true at any time, we only know the value of the ZCB exactly at a single time at the moment – the expiry date t, at which it is worth exactly $1. Plugging in these values we have

    \[Z(0,t) = \$1 \cdot {\mathbb E}\Big[\ e^{-\int_{0}^t r(t') dt' }\ |\ {\cal F}_0\ \Big]\]

So the ZCB is given by the expectation of the integral of the rate over a period of time. Since the rate is itself gaussian, and an integral is the limit of a sum, it’s not surprising that this quantity is also gaussian (it’s effectively the sum of many correlated gaussians, which is also gaussian, as discussed before), but it’s rather tricky to calculate, I’ve included the derivation at the bottom of the post to same space here. The mean and variance are given by

    \[\mathbb{E} \big[\int_s^T r(t)\ dt \ |\ {\cal F}_{\rm s}\ \big] = \big(r(s) - \theta \big) \cdot B(s,T)\ +\ \theta\cdot\big( T- s \big)\]

    \[\mathbb{V} \big[\int^T_s r(t)\ dt \ |\ {\cal F}_{\rm s}\ \big] = {\sigma^2 \over k^2} \Big( (T-s) -B(s,T) - {k\over 2} B(s,T)^2 \Big)\]

where

    \[B(s,T) = {1\over k}\big( 1 - e^{-k(T-s)}\big)\]

The ZCB is given by the expectation of the exponential of a gaussian variable – and we’ve seen on several occasions that

    \[\mathbb{E} \big[\ \exp{(x)}\ |\ x\sim {\mathbb N}(\mu, \sigma^2)\ \big] = \exp{ \Big(\mu + {1\over 2}\sigma^2 \Big) }\]

So the ZCB prices are

    \[Z(s,T) = A(s,T)\ e^{-B(s,T)\ r(s)}\]

with B(s,T) as defined above and

    \[A(s,T) = \exp{ \Big[ \Big( \theta - {\sigma^2 \over 2 k^2}\Big)\Big( B(s,T) - T + s \Big) - {\sigma^2 \over 4 k^2} B(s,T)^2 \Big]}\]

and as these expressions only depend on the rates at the initial time, we can calculate ZCB bond prices and hence forward rates for any future expiry dates at any given time if we have the current instantaneous rate.

Although we can calculate a discount curve for a given set of parameters, the Vasicek model can’t calibrate to an initial ZCB curve taken from the market, which is a serious disadvantage. There are more advanced generalisations which can, and I’ll discuss some soon, but they will use all of the same tricks and algebra that I’ve covered here.

I’ve written enough for one day here – in later posts I’ll discuss changing to the t-forward measure, in which the ZCB Z(0,t) forms the numeraire instead of a unit of currency, which simplifies many calculations, and I’ll use it to price caplets under stochastic rates, and show that these are equivalent to european options on ZCBs.

An alternate approach to the short-rate model approach discussed today which is very popular these days is the Libor Market Model (LMM) approach, in which instead of simulating the short rate and calculating the required forwards, the different forwards required are instead computed directly and in tandem – I’ll look further at this approach in another post.

 

Here is the calculation of the distribution of the integral of the instantaneous rate over the period s to T:

    \[r(t)\ = e^{-k(t-s)}\cdot r(s)\ + \theta\Big(\ 1 - e^{-k(t-s)}\ \Big) + \sigma \int_{s}^{t} e^{-k(t-u)}\cdot dW(u)\]

    \[\int_s^T r(t)\ dt = \int_s^T \Big[\ e^{-k(t-s)} r(s)\ + \theta\Big(\ 1 - e^{-k(t-s)}\ \Big) + \sigma \int_{s}^{t} e^{-k(t-u)} dW(u)\Big] dt\]

and splitting this into the terms that contribute to the expectation and the variance we have

    \[\mathbb{E} \big[\int_s^T r(t)\ dt \ |\ {\cal F}_{\rm s}\ \big] = \int_s^T \Big[\ e^{-k(t-s)} r(s)\ + \theta\Big(\ 1 - e^{-k(t-s)}\ \Big)\Big] dt\]

    \[= r(s) \int_s^T e^{-k(t-s)} \ dt + \theta \int_s^T \Big(\ 1 - e^{-k(t-s)}\ \Big) dt\]

    \[= r(s) \Big[ {1 \over -k}\ e^{-k(t-s)} \Big]^T_s + \theta \Big[ t - {1 \over -k}\ e^{-k(t-s)}\ \Big]^T_s\]

    \[= \big(r(s) - \theta \big)B(s,T)\ +\ \theta\cdot\big( T- s \big)\]

to calculate the variance, we first need to deal with the following term

    \[\int_s^T \Big[\ \sigma \int_{s}^{t} e^{-k(t-u)} dW(u)\Big] dt\]

we use stochastic integration by parts

    \[\sigma \int_s^T \int_{s}^{t} e^{-k(t-u)} dW(u)\ dt = \int_s^T e^{-kt}\big(\int_{s}^{t} e^{ku} dW(u)\big)\ dt\]

    \[= \sigma\int_s^T \big(\int_{s}^{t} e^{ku} dW(u)\big)\ d_t \Big( \int^t_s e^{-kv}dv \Big)\]

    \[= \sigma\Big[ \int_{s}^{t} e^{ku} dW(u)\cdot \int^t_s e^{-kv}dv\ \Big]^T_s - \sigma\int_s^T dW(t)\Big( e^{kt} \int^t_s e^{-kv}dv\ \Big)\]

    \[= \sigma\int_{s}^{T} e^{ku} dW(u)\cdot \int^T_s e^{-kv}dv - \sigma\int_s^T dW(t) e^{kt} \Big( \int^T_s e^{-kv}dv - \int^T_t e^{-kv}dv \Big)\]

    \[= \sigma\int_s^T e^{kt} \Big(\int^T_t e^{-kv}dv\Big) dW(t)\]

    \[= {\sigma \over k}\int_s^T \big(1 - e^{-k(T-t)}\big)\ dW(t)\]

and we’re now in a position to try and find the variance of the integral

    \[\mathbb{V} \big[\int^T_s r(t)\ dt \ |\ {\cal F}_{\rm s}\ \big] = {\mathbb E}\Big[\ \Big({\sigma \over k}\int_s^T \big(1 - e^{-k(T-t)}\big)\ dW(t)\Big)^2\ \Big]\]

    \[= {\mathbb E}\Big[\ {\sigma^2 \over k^2}\int_s^T \big(1 - e^{-k(T-t)}\big)^2\ dt\ \Big]\]

where Ito’s isometry has been used again, and several more lines of routine algebra leads to the result

    \[= {\sigma^2 \over k^2} \Big( (T-s) -B(s,T) - {k\over 2} B(s,T)^2 \Big)\]

Bootstrapping the Discount Curve from Swap Rates

Today’s post will be a short one about calculation of discount curves from swap rates. I’ve discussed both swaps and discount curves in previous posts, you should read those before this one or it might not make much sense!

Although bonds can be used to calculate discount bond prices, typically swaps are the most liquid products on the market and will go to the longest expiry times (often 80+ years for major currencies), so these are used to calculate many of the points on the discount curve [and often both of these can be done simultaneously to give better reliability].

In the previous post on swaps, I calculated the swap rate X that makes swaps zero-valued at the current time

    \[X = {Z(0,t_0) - Z(0,t_N) \over \sum_{n=0}^{N-1} \tau \cdot Z(0,t_{n+1})}\]

where the t‘s here represent the fixing dates of the swap (although payment is made at the beginning of the following period, so the n‘th period is received at t_{n+1}.

Consider the sequence of times \{t_0 = 0, t_1, t_2, \dots , t_n\} for which a sequence of swaps are quoted on the markets, with swap rates X_i for the swap running from t_0 up to t_i. We can back out the discount factor at each time as follows:

    \[X_0 = {Z(0,0) - Z(0,t_1) \over \tau \cdot Z(0,t_1)}\]

    \[\Rightarrow Z(0,t_1) = {1 \over 1 + \tau X_0 }\]

    \[X_1 = {Z(0,0) - Z(0,t_2) \over \tau \big(Z(0,t_1) + Z(0,t_2)\big) }\]

    \[\Rightarrow Z(0,t_2) = {1 - \tau \cdot X_1 \cdot Z(0,t_1) \over 1 + \tau X_1}\]

and we can see from this the general procedure, calculating another ZCB from each successive swap rate using the expression

    \[Z(0,t_i) = {1 - \sum_{i=1}^{i-1} \big(\tau \cdot X_i \cdot Z(0,t_i) \big) \over 1 + \tau X_{i-1}}\]

These swaps and ZCBs are called co-initial because they both started at the same time t_0.

Now imagine that instead the swaps X_i have the first fixing at time t_i and their final fixing at time t_n for 0 \leq i < n – such swaps are called co-terminal swaps as they start at different times but finish at the same one. Once again we can calculate the discount factors up to a constant factor, this time by working backwards:

    \[X_{n-1} = {Z(0,t_{n-1}) - Z(0,t_n) \over \tau \cdot Z(0,t_n)}\]

    \[\Rightarrow Z(0,t_{n-1}) = Z(0,t_n)\cdot\big( { 1 + \tau X_{n-1} }\big)\]

    \[X_{n-2} = {Z(0,t_{n-2}) - Z(0,t_n) \over \tau \big(Z(0,t_{n-1}) + Z(0,t_n)\big) }\]

    \[\Rightarrow Z(0,t_{n-2}) = Z(0,t_{n})\big( 1 + \tau X_{n-2} ( 1 + \tau X_{n-1}) \big)\]

and so on, the dcfs can be backed out.

To specify the exact values of co-terminal swaps, we need to know at least one dcf exactly. In general the co-initial case will also require this – I implicitly assumed that they started fixing at t=0 where we know Z(0,0) = 1, but for general co-initial swaps we would also have this issue.

Solving the Heat Equation

This post follows on from my earlier post on the BS Equation from Delta Hedging. In that post, I showed how delta hedging arguments within the BS model lead to the heat equation

{1\over 2}\sigma^2 {\partial^2 D\over \partial y^2} - {\partial D \over \partial \tau} = 0

where \inline D = C e^{-r \tau}, C(S,t) is the derivative price, \inline y = \ln S + \big( r - {1\over 2}\sigma^2 \big)\tau, and \inline \tau = T - t. In this post I’m going to solve this equation to give the BS option price for a vanilla call. S can only take positive values, but due the the logarithm y runs from \inline -\infty to \inline \infty.

In physics, this equation is usually encountered with periodic boundary conditions coming from the boundaries of the region of interest, but here we don’t have that luxury so will need to be more careful. In this post I’m going to solve the equation using the method of separation of variables, in a later post I’ll re-solve it using the method of Fourier transforms.

The heat equation is linear, which means that if \inline f(y,\tau) and \inline g(y,\tau) are both solutions, then so is the sum \inline \big( f(y,\tau) + g(y,\tau) \big). This means that if we can find a selection of possible solutions, we can combine them to match the boundary condition given by the known form of \inline D(y,\tau = 0) at expiry (coming from \inline C(S,\tau = 0) = \big( S - K \big)^+). We look for separable solutions, which are those of the form \inline f(y,\tau) = Y(y)T(\tau), in which case the derivatives in the heat equation only act on their respective terms:

{1\over 2}\sigma^2 T(\tau){\partial^2 Y(y)\over \partial y^2} - Y(y){\partial T(\tau) \over \partial \tau} = 0

{1\over 2}\sigma^2 Y''T = Y \dot{T}

{1\over 2}\sigma^2 {Y''\over Y} = {\dot{T} \over T} = \lambda^2

where a dash denotes differentiation by y and a dot differentiation by tau, In the final equation, both sides depend only on independent variables y and tau respectively – the only way this can be true is if both sides are equal to some constant, which we call \inline \lambda^2.

In general, we should consider the possibilities that \inline \lambda^2 is zero or negative, but it turns out that these won’t matter in this example so I ignore them for brevity. Two equations can be extracted, linked by \inline \lambda^2

Y'' - 2 {\lambda^2 \over \sigma^2} Y = 0

\dot{T} - \lambda^2 T = 0

Solving these equations,

Y_\lambda(y) = Y_{0,\lambda} e^{\sqrt{2}{\lambda \over \sigma }y}

T_\lambda (\tau) = e^{\lambda^2 \tau}

with any positive value of lambda allowed. We first try to find the fundamental solution of the heat equation, which is the solution that satisfies \inline f(y,\tau = 0) = \delta(y-a) where \inline \delta is the Dirac delta function, and I’ll use this to find the particular solution for a vanilla call at the end. At \inline \tau = 0 all of the functions T are 1, so these can be neglected, and we simply sum up the right number of Y terms, each with a different value of \inline \lambda which we index with n, to produce the delta function:

\sum_n Y_{0,\lambda_n} e^{\sqrt{2}{\lambda_n \over \sigma }y} = \delta(y-a)because \inline \lambda can take any value on the continuous interval \inline -\infty < \lambda < \infty, we change the sum to an integral, and also express the delta function using its integral representation and compare co-efficient to equate the two expressions

\int_{-\infty}^{\infty} Y_{0,\lambda_n} e^{\sqrt{2}{\lambda_n \over \sigma }y}dn = {1 \over 2\pi} \int_{-\infty}^{\infty} e^{in(y-a)}dn

so

Y_{0,n} = {1 \over 2\pi}e^{-in a} \quad ; \quad \lambda_n = {1 \over \sqrt{2}}\sigma i n

And plugging these into our earlier expressions for Y and T, and combining these in the same way for \inline \tau > 0 to provide the temporal behavior, we arrive at

f(y,\tau) = \int_{-\infty}^{\infty} Y_{\lambda_n}(y) T_{\lambda_n}(\tau) dn

= {1\over 2\pi} \int_{-\infty}^{\infty} e^{in(y-a)} e^{-{\sigma^2 n^2 \tau \over 2}}dn

and once again, this is a gaussian integral that we’ll tackle by completing the square…

= {1\over 2\pi} \int_{-\infty}^{\infty} \exp\Big\{ -{1 \over 2} \sigma^2 \tau \Big( n^2 - 2 {i(y-a)\over \sigma^2 \tau} + \big( {i(y-a)\over \sigma^2 \tau}\big)^2-\big( {i(y-a)\over \sigma^2 \tau}\big)^2 \Big) \Big\} dn

= {1\over 2\pi} \int_{-\infty}^{\infty} \exp\Big\{ -{1 \over 2} \sigma^2 \tau \Big(n- \big( {i(y-a)\over \sigma}\big) \Big)^2 \Big\} \cdot \exp{-(y-a)^2 \over 2 \sigma^2 \tau} dn

= {1\over 2\pi}\cdot\sqrt{2\pi \over \sigma^2 \tau}\cdot \exp{-(y-a)^2 \over 2 \sigma^2 \tau}

= {1\over \sqrt{2\pi \sigma^2 \tau}}\cdot \exp{-(y-a)^2 \over 2 \sigma^2 \tau}

and this is the fundamental solution given at the end of the first post and also at wikipedia. Note that although we used separation of variables to find the forms of the functions Y and T, the separability hasn’t been carried through the summation of many different T-Y combinations – the fundamental solution that we’ve arrived at ISN’T separable!

The fundamental solution is the time and space evolution of a function that is initially a delta-function, such that the evolution always obeys the heat equation. We can use the following property of the delta function to construct specific solutions for particular payoff functions:

f(y) = \int_{-\infty}^{\infty} \delta(y-a) f(a) da

So, because the fundamental solution \inline f(y,\tau=0) = \delta(y-a), we can construct the payoff function using this integral at \inline \tau = 0

D(y,\tau=0) = \int_{-\infty}^{\infty} f(a,\tau=0) D(a,\tau=0) da

and consequently, at later times

D(y,\tau) = \int_{-\infty}^{\infty} f(a,\tau) D(a,\tau=0) daand we can solve this for the specific case that \inline C(S,\tau = 0) = \big( S(\tau=0) - K \big)^+ as follows (denoting a as y’)

D(y,\tau) = {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{-\infty}^{\infty} \exp\Big\{ - {(y'-y)^2 \over 2\sigma^2 \tau} \Big\} \cdot \big( e^{y'} - K \big)^+ dy'

The second factor is only non-zero when \inline y' > \ln K = y_L, so we can absorb this condition into the lower limit

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {(y'-y)^2 \over 2\sigma^2 \tau} \Big\} \cdot \big( e^{y'} - K \big) dy'and splitting this into two integrals, which we do in turn

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} K\exp\Big\{ - {(y'-y)^2 \over 2\sigma^2 \tau} \Big\} dy' = e^{-r\tau}\cdot K\cdot \Phi\Big({y - y_L \over \sigma \sqrt{\tau}}\Big)

and

{e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {(y'-y)^2 \over 2\sigma^2 \tau} \Big\} \cdot e^{y'} dy'

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {1\over 2\sigma^2 \tau}\Big( y'^2 + y^2 - 2y'y -2y'\sigma^2\tau \pm ( y + \sigma^2\tau )^2 \Big) \Big\} dy'

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {1\over 2\sigma^2 \tau}\Big( \big(y'-(y+\sigma^2 \tau)\big)^2 - \big(y + \sigma^2 \tau\big)^2 + y^2\Big) \Big\} dy'

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {1\over 2\sigma^2 \tau}\Big( \big(y'-(y+\sigma^2 \tau)\big)^2 - \sigma^2\tau\big(2y + \sigma^2 \tau\big)\Big) \Big\} dy'

= {e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {\big(y'-(y+\sigma^2 \tau)\big)^2\over 2\sigma^2 \tau}\Big\}\cdot S e^{r\tau} dy'

= e^{-r\tau} \cdot F \cdot \Phi \Big( {y - y_L + \sigma^2 \tau \over \sigma \sqrt{\tau}} \Big)

where as usual, \inline \Phi(x) is the cumulative normal distribution; and we’ve used the expression \inline y = \ln S + (r - {1\over 2}\sigma^2)\tau and also \inline y + {1 \over 2}\sigma^2 \tau = S e^{r\tau} = F. Putting all of this together, we finally arrive at the Black-Scholes expression for a vanilla call price,

{e^{-r\tau} \over \sqrt{2\pi\sigma^2\tau}}\int_{y_L}^{\infty} \exp\Big\{ - {(y'-y)^2 \over 2\sigma^2 \tau} \Big\} \cdot \big( e^{y'} - K \big) dy'

= e^{-r\tau} \Big( F \cdot \Phi(d_1) - K\cdot \Phi(d_2)\Big)

 

 

Wow! Finally! Hopefully this will convince you that the method of risk-neutral valuation is a bit more straight-forward!

Digital Options

Today I’m going to talk about the valuation of another type of option, the digital (or binary) option. This can be seen as a bit of a case study, as I’ll present the option payoff and the analytical price and greeks under BS assumptions, and give add-ons to allow pricing with the MONTE CARLO pricer. I’ve also updated the ANALYTICAL pricer to calculate the price and greeks of these options.

Digital options are very straight-forward, they are written on an underlying S and expire at a particular date t, at which point digital calls pay $1 if S is greater than a certain strike K or $0 if it is below that, and digital puts pay the reverse – ie. the payoff is

P_{\rm dig\ call} = \Big\{ \ \begin{matrix} \$1 \quad S \geq K\\ \$0 \quad S < K\end{matrix}

We can calculate the price exactly in the BS approximation using the same method that I used to calculate vanilla option prices by risk-neutral valuation as follows

C_{\rm dig\ call}(0) = \delta(0,t)\ {\mathbb E}[ C_{\rm dig\ call}(t) ]

where \inline C_{\rm dig\ call}(t') is the price of the option at time \inline t', and we know that this must converge to the payoff as \inline t \to t', so \inline C_{\rm dig\ call}(t) = P_{\rm dig\ call}

= \delta(0,t)\ \int^{\infty}_{-\infty} P_{\rm dig\ call} \cdot e^{-{1\over 2}x^2} dx

\inline P_{\rm dig\ call} is zero if \inline S = S_0 e^{(r - {1\over 2}\sigma^2)t + \sigma \sqrt{t} x} < K which corresponds to \inline x < {\ln{K \over S_0} - (r-{1\over 2}\sigma^2)t \over \sigma \sqrt{t}} = -d_2

= \delta(0,t)\ \cdot \$1 \cdot \int^{\infty}_{-d_2} e^{-{1\over 2}x^2} dx

= \$ \delta(0,t)\cdot \Phi(d_2)

where

d_1 = {\ln{S\over K}+(r+{1\over 2}\sigma^2)t \over \sigma \sqrt{t}} \quad ;\quad d_2 = {\ln{S\over K}+(r-{1\over 2}\sigma^2)t \over \sigma \sqrt{t}}

Since we have an analytical price, we can also calculate an expression for the GREEKS of this option by differentiating by the various parameters that appear in the price. Analytical expressions for a digital call’s greeks are:

\Delta = {\partial C \over \partial S} = e^{-rt}\cdot {\phi(d_2)\over S\sigma \sqrt{t}}

\nu = {\partial C \over \partial \sigma} = -e^{-rt}\cdot {d_1 \ \phi(d_2)\over \sigma}

\gamma = {\partial^2 C\over \partial S^2} = -e^{-rt}\cdot {d_1 \ \phi(d_2)\over S^2 \sigma^2 t} = -e^{-rt}\cdot {d_1\ \phi(d_1)\over S K \sigma^2 t}

{\rm Vanna} = {\partial^2 C \over \partial S \partial \sigma} =e^{-rt}\cdot {\phi(d_2)\over S \sigma^2 \sqrt{t}}\Big( d_1 d_2 - 1 \Big)

{\rm Volga} = {\partial^2 C \over \partial \sigma^2} = e^{-rt}\cdot {\phi(d_2)\over \sigma^2}\cdot \Big( d_1 + d_2 - d_1^2 d_2 \Big)

where

\phi(d_1) = {1\over \sqrt{2 \pi}}\ {e^{-{1\over 2}d_1^2} }

Holding a binary put and a binary call with the same strike is just the same as holding a zero-coupon bond, since we are guaranteed to receive $1 wherever the spot ends up, so the price of a binary put must be

e^{-rt} = e^{-rt}\cdot \Phi(d_2) + C_{\rm dig\ put}(t'=0)\quad \quad \quad \quad [1]

C_{\rm dig\ put}(0) = e^{-rt} \cdot \Phi(-d_2)

Moreover, differentiating equation [1] above shows that the greeks of a digital put are simply the negative of the greeks of a digital call with the same strike.

Graphs of these are shown for a typical binary option in the following graphs. Unlike vanilla options, these option prices aren’t monotonic in volatility: if they’re in-the-money, increasing vol will actually DECREASE the price since it makes them more likely to end out-of-the-money!

Price and first-order greeks for a digital call option.
Price and first-order greeks for a digital call option.
DigitalGreeks2
Second-order greeks for a digital call option. Greeks for digital puts are simply the negative of these values

One final point on pricing, note that the payoff of a digital call is the negative of the derivative of a vanilla call payoff wrt. strike; and the payoff of a digital put is the positive of the derivative of a vanilla put payoff wrt. strike. This means that any binary greek can be calculated from the corresponding vanilla greek as follows

\Omega_{\rm dig\ call} = -{\partial \Omega_{\rm call}\over \partial K}

\Omega_{\rm dig\ put} = {\partial \Omega_{\rm put}\over \partial K}

where here \inline \Omega represents a general greek.

If you haven’t yet installed the MONTE CARLO pricer, you can find some instructions for doing so in a previous post. The following links give the header and source files for binary calls and puts which can be dropped in to the project in your C++ development environment

These will register the option types with the option factory and allow monte carlo pricing of the options (so far, all of the options in the factory also have analytical expressions, but I’ll soon present some options that can only be priced by Monte Carlo).

Interview Quesions III

Today’s question will test some of the statistics and correlation I’ve discussed in the last couple of months. Assume throughout that x\sim {\mathbb N}(0,1) and y\sim {\mathbb N}(0,1) are jointly normally distributed such that {\mathbb E}[x \cdot y] = \rho

a) Calculate {\mathbb E}[\ e^x \ ]
b) Calculate {\mathbb E}[\ e^x \ | \ y = b\ ]

The first expectation is of a lognormal variate, and the second is of a lognormal variate conditional on some earlier value of the variate having been a particular value – these are very typical of the sorts of quantities that a quant deals with every day, so the solution will be quite instructive! Before reading the solution have a go at each one, the following posts may be useful: SDEs pt. 1, SDEs pt. 2, Results for Common Distributions

a) Here, we use the standard result for expectations

    \begin{align*} {\mathbb E}[ \ e^x \ ] &= \int^{\infty}_{-\infty} e^x \cdot f(x) \ dx \nonumber \\ \nonumber \\ &= {1 \over \sqrt{2 \pi}}\int^{\infty}_{-\infty} e^x \cdot e^{-{1 \over 2}x^2} \ dx \nonumber \\ \nonumber \\ &= {1 \over \sqrt{2 \pi}}\int^{\infty}_{-\infty} \exp\Bigl( -{1\over 2}\Bigl[ x^2 - 2x + 1 - 1 \Bigr] \Bigr) \ dx \nonumber \\ \nonumber \\ &= {e^{1\over 2} \over \sqrt{2 \pi}}\int^{\infty}_{-\infty} \exp\Bigl( -{1\over 2} (x-1)^2 \Bigr) \ dx \nonumber \\ \nonumber \\ &= e^{1\over 2} \nonumber \end{align}

b) This one is a little tougher, so first of all I’ll discuss what it means and some possible plans of attack. We want to calculate the expectation of e^x, given that y takes a value of b. Of course, if x and y were independent, this wouldn’t make any difference and the result would be the same. However, because they are correlated, the realised value of y will have an effect on the distribution of x.

To demonstrate this, I’ve plotted a few scatter-graphs illustrating the effect of specifying y on x, with x and y uncorrelated and then becoming increasing more correlated.

When x and y are uncorrelated, the realised value of y doesn't affect the distribution for x, which is still normally distributed around zero
When x and y are uncorrelated, the realised value of y doesn’t affect the distribution for x, which is still normally distributed around zero
When x and y are correlated, the realised value of y has an effect on the distribution of x, which is no longer centered on zero and has a smaller variance
When x and y are correlated, the realised value of y has an effect on the distribution of x, which is no longer centered on zero and has a smaller variance
When the correlation of x and y becomes high, the value of x is almost completely determined by y. Now, if y is specified then x is tightly centered around a value far from zero
When the correlation of x and y becomes high, the value of x is almost completely determined by y. Now, if y is specified then x is tightly centered around a value far from zero

The simplest way of attempting this calculation is to use the result for jointly normal variates given in an earlier post, which says that if x and y have correlation \rho, we can express x in terms of y and a new variate z \sim {\mathbb N}(0,1) which is uncorrelated with y

    \[x = \rho\cdot y + \sqrt{1 - \rho^2}\cdot z\]

so

    \[(\ e^x\ |\ y = b\ ) = e^{\rho y + \sqrt{1-\rho^2}z} = e^{\rho b}\cdot e^{\sqrt{1-\rho^2}z}\]

Since the value of y is already determined (ie. y = b), I’ve separated this term out and the only thing I have to calculate is the expectation of the second term in z. Since y and z are independent, we can calculate the expectation of the z, which is the same process as before but featuring slightly more complicated pre-factors

    \begin{align*} {\mathbb E}[ \ e^x \ |\ y = b\ ] &= e^{\rho b} \int^{\infty}_{-\infty} e^{\sqrt{1-\rho^2}z} \cdot f(z) \ dz \nonumber \\ \nonumber \\ &= {e^{\rho b} \over \sqrt{2 \pi}}\int^{\infty}_{-\infty} e^{\sqrt{1-\rho^2}z} \cdot e^{-{1 \over 2}z^2} \ dz \nonumber \\ \nonumber \\ &= {e^{\rho b} \over \sqrt{2 \pi}}\int^{\infty}_{-\infty} \exp\Bigl( -{1\over 2}\Bigl[ z^2 - 2\sqrt{1-\rho^2}z \nonumber \\ & \quad \quad \quad \quad \quad \quad \quad \quad + (1-\rho^2) - (1-\rho^2) \Bigr] \Bigr) \ dz \nonumber \\ \nonumber \\ &= {e^{\rho b} \over \sqrt{2 \pi}} e^{{1 \over 2} (1-\rho^2)} \int^{\infty}_{-\infty} \exp\Bigl( -{1\over 2} \bigl(z-\sqrt{1-\rho^2}\bigr)^2 \Bigr) \ dz \nonumber \\ \nonumber \\ &= e^{{1\over 2}(1-\rho^2)+\rho b} \nonumber \end{align}

We can check the limiting values of this – if \rho = 0 then x and y are independent [this is not a general result by the way – see wikipedia for example – but it IS true for jointly normally distributed variables], in this case {\mathbb E}[ \ e^x \ |\ y = b\ ] = e^{0.5} just as above. If \rho = \pm 1, {\mathbb E} [\ e^x \ |\ y=b\ ] = e^{\pm b}, which also makes sense since in this case x = \pm y = \pm b, so y fully determines the expectation of x.

The more general way to solve this is to use the full 2D joint normal distribution as given in the previous post mentioned before,

    \[f(x,y) = {1 \over {2\pi \sqrt{1-\rho^2}}} \exp{\Bigl(-{1 \over 2(1-\rho^2)}(x^2 + y^2 - 2\rho xy)\Bigr)}\]

This is the joint probability function of x and y, but it’s not quite what we need – the expectation we are trying to calculate is

    \[{\mathbb E}[ \ e^x \ |\ y \ ] = \int^{\infty}_{-\infty} e^x \cdot f(\ x \ |\ y \ ) \ dx\]

So we need to calculate the conditional expectation of x given y, for which we need Bayes’ theorem

    \[f(x,y) = f( \ x \ | \ y \ ) \cdot f(y)\]

Putting this together, we have

    \[{\mathbb E}[ \ e^x \ |\ y = b\ ] = \int^{\infty}_{-\infty} e^x \cdot {f(x,y) \over f(y)}\ dx\]

    \[={1 \over {2\pi \sqrt{1-\rho^2}}} \int^{\infty}_{-\infty} { e^x \over e^{-{1\over 2}y^2}} \exp{\Bigl(-{1 \over 2(1-\rho^2)}(x^2 + y^2 - 2\rho xy)\Bigr)} dx\]

    \[={e^{{1\over 2}b^2} \over {2\pi \sqrt{1-\rho^2}}} \int^{\infty}_{-\infty} e^x \exp{\Bigl(-{1 \over 2(1-\rho^2)}(x^2 + b^2 - 2\rho xb)\Bigr)} dx\]

This integral is left as an exercise to the reader, but it is very similar to those given above and should give the same answer as the previous expression for {\mathbb E}[\ e^x \ | \ y = b\ ] after some simplification!