Black-Scholes Call and Put Equations and Comparative Statics

Latest Version: April 23, 2019
This Notebook describes how to use Mathematica to perform generate graphs of the so-called option "Greeks". Suggestions concerning ways to improve this notebook, both in terms of its programming or content, are quite welcome and should be sent to James_Garven@baylor.edu.

Summary of Black-Scholes Option Pricing Formulas and Comparative Statics

The Black-Scholes formula for the value of a European call option is:
c=SN((Subscript[d, 1]))−K(e^(−rτ))N((Subscript[d, 2])), where (Subscript[d, 1])==((ln(S/K)+(r+.5(σ^2))τ)/(σ(Sqrt[τ]))),  (Subscript[d, 2])=(Subscript[d, 1])−σ(Sqrt[τ]), and τ=T−t; i.e., τ represents the time to expiration for the option.
The Black-Scholes formula for the value of a European put option is:
p=K(e^(−rτ))N(-(Subscript[d, 2])) - SN(-(Subscript[d, 1])) .
Thus, call and put option prices are functions of 1) the current asset price S, 2) the exercise price K, 3) the rate of interest r, 4) the time to expiration t, and 5) volatility σ (option prices are functions of the dividend yield, but in the numerical simulations which follow, we’ll assume that the underlying asset does not pay dividends).  
Thus, the comparative statics of call and put option prices with respect to S, K, r, t, and σ are as follows:
1. Asset Price (“Delta”): ∂c/∂S=N((Subscript[d, 1]))>0 and ∂p/∂S=−N(−(Subscript[d, 1]))<0.
2. Exercise Price: ∂c/∂K=−(e^(−rτ))N((Subscript[d, 2]))<0 and  ∂p/∂K=(e^(−rτ))N(−(Subscript[d, 2]))>0.
3. Interest Rate (“rho”): ∂c/∂r=τK(e^(−rτ))N((Subscript[d, 2]))>0 and ∂p/∂r=−τK(e^(−rτ))N(−(Subscript[d, 2]))<0.
4. Time to Expiration (“theta”): ∂c/∂t=−Sn((Subscript[d, 1]))((.5σ)/(Sqrt[τ]))−rK(e^(−rτ))N((Subscript[d, 2]))<0 and ∂p/∂t=−Sn((Subscript[d, 1]))((.5σ)/(Sqrt[τ]))+rK(e^(−rτ))N(−(Subscript[d, 2]))<>0.
5. Volatility (“vega”): ∂c/∂σ=Sn((Subscript[d, 1]))(Sqrt[τ]); ∂p/∂σ=Sn(−(Subscript[d, 1]))(Sqrt[τ]). Note that ∂c/∂σ=∂p/∂σ>0.
6. Gamma: (∂^2)c/∂(S^2)=(∂^2)p/∂(S^2)=n((Subscript[d, 1]))/Sσ(Sqrt[τ])>0.

Numerical Simulations

We begin our numerical analysis by coding the inputs for the Black-Scholes call and put option equations and comparative statics:

d1[S_, K_, r_, t_, q_, σ_] := ((r - q)*t + Log[S/K])/(σ*Sqrt[t]) + (σ*Sqrt[t])/2; 
d2[S_, K_, r_, t_, q_, σ_] := d1[S, K, r, t, q, σ] - σ*Sqrt[t]; 
minusd1[S_, K_, r_, t_, q_, σ_] := d1[S, K, r, t, q, σ]*(-1); 
minusd2[S_, K_, r_, t_, q_, σ_] := d2[S, K, r, t, q, σ]*(-1);
Normsdist[K_] := 1/2*(1 + Erf[K/Sqrt[2]]); 
BlackScholesCall[S_, K_, r_, t_, q_, σ_] := S*Exp[-q*t]*Normsdist[d1[S, K, r, t, q, σ]] - K*Exp[-r*t]*Normsdist[d2[S, K, r, t, q, σ]]; 
BlackScholesPut[S_, K_, r_, t_, q_, σ_] := K*Exp[-r*t]*Normsdist[minusd2[S, K, r, t, q, σ]] - S*Exp[-q*t]*Normsdist[minusd1[S, K, r, t, q, σ]]; 
dcdk[S_, K_, r_, t_, q_, σ_] := -Exp[-r*t]*Normsdist[d2[S, K, r, t, q, σ]];
dputdk[S_, K_, r_, t_, q_, σ_] := Exp[-r*t]*Normsdist[minusd2[S, K, r, t, q, σ]];
dcalldr[S_, K_, r_, t_, q_, σ_] := t*K*Exp[-r*t]Normsdist[d2[S, K, r, t, q, σ]];
dputdr[S_, K_, r_, t_, q_, σ_] := -t*K*Exp[-r*t]Normsdist[minusd2[S, K, r, t, q, σ]];
nd1[S_, K_, r_, t_, q_, σ_] := Exp[-.5*d1[S, K, r, t, q, σ]^2]/Sqrt[2*π];
nd1minus[S_, K_, r_, t_, q_, σ_] := Exp[-.5*minusd1[S, K, r, t, q, σ]^2]/Sqrt[2*π];
nd2[S_, K_, r_, t_, q_, σ_] := Exp[-.5*d2[S, K, r, t, q, σ]^2]/Sqrt[2*π];
gamma[S_, K_, r_, t_, q_, σ_] := nd1[S, K, r, t, q, σ]/S*σ/Sqrt[t];
exercisegamma[S_, K_, r_, t_, q_, σ_] := Exp[-r*t]nd2[S, K, r, t, q, σ]/S*σ/Sqrt[t];
calltheta[S_, K_, r_, t_, q_, σ_] := -r*K*Exp[-r*t]Normsdist[d2[S, K, r, t, q, σ]]-S*nd1[S, K, r, t, q, σ]*.5*σ/Sqrt[t];
puttheta[S_, K_, r_, t_, q_, σ_] := r*K*Exp[-r*t]Normsdist[minusd2[S, K, r, t, q, σ]]-S*nd1[S, K, r, t, q, σ]*.5*σ/Sqrt[t];
vega[S_, K_, r_, t_, q_, σ_] := S*nd1[S, K, r, t, q, σ]*Sqrt[t];
putvega[S_, K_, r_, t_, q_, σ_] := S*nd1minus[S, K, r, t, q, σ]*Sqrt[t];

“Base Case” values for the call and put option

Thus S_ corresponds to the current value of the underlying, K_ corresponds to the strike price, r_ corresponds to the interest rate, t_ corresponds to the time to expiration, q_ corresponds to the dividend rate, and σ_ corresponds to volatility. The base case values for the option parameters will be S_ = 50, K_ = 50, r_ = 5%, t_ = 1 year, q_ = 0, and σ_=30% (these are the parameter values used in Hull’s chapter entitled “The Black-Scholes-Merton Model”).  Thus, the “base case” values for the call and put option are as follows:

BlackScholesCall[50, 50, .05, 1, 0, .3]

BlackScholesPut[50, 50, .05, 1, 0, .3]

Normsdist[d1[50,50,0.05,1,0,0.3]]

-Normsdist[-d1[50,50,0.05,1,0,0.3]]

gamma[50,50,0.05,1,0,0.3]

-Exp[-.05(1)]*Normsdist[d2[50,50,0.05,1,0,0.3]]

Exp[-.05(1)]*Normsdist[-d2[50,50,0.05,1,0,0.3]]

1* 50*Exp[-.05(1)]*Normsdist[d2[50,50,0.05,1,0,0.3]]

-1* 50*Exp[-.05(1)]*Normsdist[-d2[50,50,0.05,1,0,0.3]]

calltheta[50,50,0.05,1,0,0.3]

puttheta[50,50,0.05,1,0,0.3]

vega[50,50,0.05,1,0,0.3]

putvega[50,50,0.05,1,0,0.3]

Relationship between Option Deltas and the price of the underlying

First, we test the relationship between the call option delta and the price of the underlying:

Plot[Normsdist[d1[underlying,50,0.05`,1,0,0.3`]],{underlying,0,100},PlotLabel->"Call Delta as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Delta"}]

As one would intuitively expect, as the price of the underlying increases, then the hedge ratio converges toward 1.  This makes sense because if the option is in-the-money, then the ratio of changes in the option value with respect to changes in the value of the underlying become highly positively correlated.  On the other hand, as the price of the underlying diminishes, then the value of the option becomes much less correlated with the value of the underlying; hence the (near zero) hedge ratio in such cases.
One can intuitively infer the "Gamma" for the call option visually by noting that Gamma = (((∂^2)c)/(∂(S^2))); i.e., it measures how the hedge ratio changes as the price of the underlying varies. Gamma has a very small value for options which are either deeply out-of-the-money or deeply in-the-money (note that the slope is quite flat in these cases), but for options that are near the money, Gamma is highly sensitive to changes in the price of the underlying asset.  
Next, we test the relationship between the put hedge ratio and the price of the underlying:

Plot[-Normsdist[-d1[underlying,50,0.05`,1,0,0.3`]],{underlying,0,100},PlotLabel->"Put Delta as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Delta"}]

Not surprisingly, the graph for the put looks very much like the graph for the call, except the y axis values are negative. Of course, we expect this since delta for the put option (-N(-(Subscript[d, 1]))) is equal to the call option's delta (N((Subscript[d, 1]))) minus 1!  Here, the hedge ratio for the put option converges toward 0 for deeply out-of-the-money puts (i.e., when the price of the underlying is high), and toward -1 for deeply  The same point made concerning Gamma for call options applies to put options; i.e., Gamma has a very small value for options which are either deeply out-of-the-money or deeply in-the-money (note that the slope is quite flat in these cases), but for options that are near the money, gamma is highly sensitive to changes in the price of the underlying asset.

Relationship between Option Gammas and the price of the underlying

Plot[gamma[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Gamma as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Gamma"}]

Here, Gamma (for both the call and put) is maximized at a stock price slightly above $40.  As noted earlier in our analysis of Delta, Gamma has a very small value for options which are either deeply out-of-the-money or deeply in-the-money, but for options that are near the money, Gamma is highly sensitive to changes in the price of the underlying asset.

Relationship between ∂c/∂K and ∂p/∂K and the price of the underlying

Plot[dcdk[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"∂c/∂K as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","∂c/∂K"}]

As in the case of the option Gammas, we find that the sensitivity of the option price with respect to changes in the exercise price is particularly high for near the money options but rather insensitive in the case of out-of-the-money options. Also, note that the net difference between  ((∂p)/(∂K))and ((∂c)/(∂K))  is equal to (e^(-rt )).

Relationship between Option Rhos and the price of the underlying

Plot[dcalldr[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Call Rho as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Rho"}]

Note that call option value is not sensitive to interest rates when the stock price is low (i.e., when the call is out-of-the-money); however, call option value is very sensitive to interest rates for high stock prices (i.e., when the call is in-the-money).

Plot[dputdr[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Put Rho as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Rho"}]

Not surprisingly, put option value is highly sensitive to interest rates when the stock price is low (i.e., when the put is in-the-money); however, put option value is not sensitive to interest rates for high stock prices.  (i.e., when the put is out-of-the-money).

Relationship between Option Thetas and the price of the underlying

Plot[calltheta[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Call Theta as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Theta"}]

Theta is always negative for a call option.  This is because, as time passes with all else remaining the same, the call option becomes less valuable.  However, the rate at which this occurs depends upon the moneyness of the option.  The variation of theta with stock price for a call option is shown above.  When the stock price is very low, theta is close to zero.  For an at-the-money call option, theta is large and negative.  As the stock price becomes larger, theta converges in value toward  -r(e^(-rt ))KN((Subscript[d, 2])) ;  Since the call theta is ∂c/∂t=−Sn((Subscript[d, 1]))((.5σ)/(Sqrt[τ]))−rK(e^(−rτ))N((Subscript[d, 2])), the first term vanishes since n((Subscript[d, 1]))->0 in such cases!

Plot[puttheta[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Put Theta as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Theta"}]

The pattern for the put theta is similar to the pattern for the call theta, only here we find that when the stock price is very low, theta is strongly positive.  This is because with the passage of time, the odds of a large move to a high price diminishes, so the first term on the right-hand side of the put theta equation (∂p/∂t=−Sn((Subscript[d, 1]))((.5σ)/(Sqrt[τ]))+rK(e^(−rτ))N(−(Subscript[d, 2]))) becomes negligible (since n((Subscript[d, 1]))->0 in such cases).  However, at higher stock prices  (in this case, above $38 or so), theta turns negative as the −Sn((Subscript[d, 1]))((.5σ)/(Sqrt[τ])) term dominates the rK(e^(−rτ))N(−(Subscript[d, 2])) term, with the difference between the two terms converging toward zero as the stock price becomes very high. 
Next, let’s examine how the rate of time decay is related to the passage of time for an at-the-money (ATM) call option:

Plot[calltheta[50,50,0.05`,t,0,0.3`],{t,.25,2},PlotLabel->"Call Theta for ATM option",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Time","Theta (ATM)"}]

The above graph shows that the call option’s theta becomes more negative with the passage of time.  Thus, the rate at which time decay occurs is faster for short lived options than it is for longer lived options.  Furthermore, the rate at which time decay changes with the passage of time is increasing, which implies that short-lived options (particularly options that are near or in the money) will lose value most rapidly as the expiration date draws close.

Next, let’s examine how the rate of time decay is related to the passage of time for the put option; first consider at-the-money (ATM) put (where the current stock price is $50):

Plot[puttheta[50,50,0.05`,t,0,0.3`],{t,.25,2},PlotLabel->"Put Theta for ATM option",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Time","Theta (ATM)"}]

The pattern for the ATM put theta is very similar to the pattern for the ATM  call theta, only somewhat less attenuated.  However, since the put theta is positive for deeply in-the-money (ITM) puts, let’s draw the same graph for an ITM put (where the current stock price is $30):

Plot[puttheta[30,50,0.05`,t,0,0.3`],{t,.25,2},PlotLabel->"Put Theta for ITM option",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Time","Theta (ITM)"}]

Here, the pattern for the ITM put theta is the reverse of the ATM put theta.  Specifically, rather than decay with the passage of time, the ITM put theta expands with the passage of time.  The reason for this is that as the time to maturity decreases, this makes it much more likely that the put will end up expiring in the money.

Relationship between Option Vegas and the price of the underlying

Plot[vega[underlying,50,0.05`,1,0,0.3`],{underlying,0,100},PlotLabel->"Option Vegas as a Function of Stock Price",Frame->True,GridLines->Automatic,PlotRange->All,FrameLabel->{"Stock Price","Vega"}]

Note that vega is maximized when either the call or the put is at-the-money (i.e., S = 50).