In [1]:
using Dates, Distributions, Roots, QuadGK
include("jlFiles/printmat.jl")
Out[1]:
In [2]:
using Plots
#pyplot(size=(600,400))
gr(size=(480,320))
default(fmt = :svg)
The Black-Scholes formula for a European call option on an asset with a continuous dividend rate $\delta$ is
$C =e^{-\delta m}S\Phi(d_{1}) - e^{-ym}K\Phi(d_{2})$, where
$d_{1} =\frac{\ln(S/K)+(y-\delta+\sigma^{2}/2)m}{\sigma\sqrt{m}} \ \text{ and } \ d_{2}=d_{1}-\sigma\sqrt{m}$
and where $\Phi(d)$ denotes the probability of $x\leq d$ when $x$ has an $N(0,1)$ distribution.
In [3]:
Φ(x)= cdf(Normal(0,1),x) #Calculate Pr(z<=x) for N(0,1) variable z
"""
Calculate Black-Scholes European call or put option price, continuous dividends of δ
"""
function OptionBlackSPs(S,K,m,y,σ,δ=0,PutIt=false)
d1 = ( log(S/K) + (y-δ+0.5*σ^2)*m ) / (σ*sqrt(m))
d2 = d1 - σ*sqrt(m)
c = exp(-δ*m)*S*Φ(d1) - K*exp(-y*m)*Φ(d2)
if PutIt
price = c - exp(-δ*m)*S + exp(-y*m)*K
else
price = c
end
return price
end
Out[3]:
In [4]:
(S,K,m,y,σ) = (42,42,0.5,0.05,0.2)
C = OptionBlackSPs(S,K,m,y,σ)
printlnPs("\n","call price at K=$K: ",C)
P = OptionBlackSPs(S,K,m,y,σ,0,true)
printlnPs("\n","put price at K=$K: ",P)
δ = 0.05
Pδ = OptionBlackSPs(S,K,m,y,σ,δ,true)
printlnPs("\n","put price at K=$K when δ=$δ: ",Pδ)
In [5]:
K_range = 30:60 #different strike prices
C_K_range = OptionBlackSPs.(S,K_range,m,y,σ)
plot( K_range,C_K_range,
linecolor = :red,
legend = false,
title = "Black-Scholes call option price",
xlabel = "strike price",
ylabel = "option price" )
Out[5]:
In [6]:
S_range = 30:60 #different spot prices
C_S_range = OptionBlackSPs.(S_range,K,m,y,σ)
plot( S_range,C_S_range,
linecolor = :blue,
legend = false,
title = "Black-Scholes call option price",
xlabel = "current asset price",
ylabel = "option price" )
Out[6]:
The Black-Scholes option price is an increasing function of the volatility ($\sigma$), as illustrated below.
In the subsequent cell, we use an observed option price and solve the BS formula for $\sigma$. This is the "implied volatility".
In [7]:
σ_range = 0.01:0.01:0.3 #S,K,m,y are scalars, σ_range is a vector
C_σ_range = OptionBlackSPs.(S,K,m,y,σ_range) #at different sigmas
plot( σ_range,C_σ_range,
linecolor = :red,
label = "BS price",
title = "Black-Scholes call option price",
xlabel = "volatility (sigma)",
ylabel = "option price",
legend = :top )
hline!( [C;C+0.5;C-0.5],linetype=:hline,linecolor=:black,line=(:dash,0.5),label="data on prices")
Out[7]:
In [8]:
println("Compare the following results with the previous graph\n")
iv = find_zero(sigma->OptionBlackSPs(S,K,m,y,sigma)-C,(-1,1))
printlnPs("Implied volatility (same as above?): ",iv)
iv_a = find_zero(sigma->OptionBlackSPs(S,K,m,y,sigma)-(C+0.5),(-1,1))
printlnPs("Implied volatility (overpriced option): ",iv_a)
iv_b = find_zero(sigma->OptionBlackSPs(S,K,m,y,sigma)-(C-0.5),(-1,1))
printlnPs("Implied volatility (underpriced option): ",iv_b)
The price of a European a call option is
$ C=e^{-ym}\text{E}^{\ast}\max(0,S_{m}-K), $
which can be written
$ C=e^{-ym}\int_{0}^{\infty}\max(0,S_{m}-K) f^{\ast}(S_{m})dS_{m}. $
where $f^{\ast}(S_{m})$ is the risk neutral density function of the asset price at expiration ($S_{m}$). (The spot price cannot go below zero.)
In the Black-Scholes model, the risk neutral distribution of $\ln S_{m}$ is
$ \ln S_{m}\sim^{\ast}\text{N}(\ln S+ym-\sigma^{2}m/2,\sigma^{2}m), $
where $S$ is the current asset price. This means that $f^{\ast}()$ is the pdf of a lognormally distributed variable with parameters $\ln S+ym-\sigma^{2}m/2$ and $\sigma^{2}m$. However, the Distributions package has another convention and wants the square root of the second parameter.
The numerical integration is done by the QuadGK package.
In [9]:
"""
BSintegrand(Sm,S,K,y,m,σ)
Constructs the integrand for the Black-Scholes call price
"""
function BSintegrand(Sm,S,K,y,m,σ)
μ = log(S) + y*m - m*σ^2/2 #"mean" of log-normal
f = pdf(LogNormal(μ,sqrt(m)*σ),Sm) #log-normal pdf
z = exp(-y*m)*max(0,Sm-K)*f
return z
end
Out[9]:
In [10]:
Sₘ = (K-10):0.25:(K+40) #possible outcomes of underlying price
txt = text("Integrate this to get the\nBlack-Scholes call option price",8,:left)
plot( Sₘ,BSintegrand.(Sₘ,S,K,y,m,σ),
linecolor = :red,
ylim = (-0.01,0.25),
legend = false,
title = "Integrand for B-S call price",
xlabel = "Asset price at expiration",
ylabel = "integrand",
annotation = (55,0.22,txt) )
Out[10]:
In [11]:
C1, = QuadGK.quadgk(x->BSintegrand(x,S,K,y,m,σ),0,Inf) #numerical integration
printlnPs("\nCall price according to integration of risk-neutral dist: ",C1)
printlnPs("\nCompare with the result from the BS formula (from above): ",C)
The next few cells calculate the option price according to binomial model with a CRRA calibration where
$ u=e^{\sigma\sqrt{h}}, d=e^{-\sigma\sqrt{h}} \text{ and } p=\frac{e^{yh}-d}{u-d}. $
This is done repeatedly, using more and more time steps ($n$) so that $h=m/n$ where $m$ is the fixed time to expiration.
In [12]:
include("jlFiles/OptionsCalculations.jl")
Out[12]:
In [13]:
#(S,K,m,y,σ) = (42,42,0.5,0.05,0.2) #these parameters were defined before
nMax = 200
CeM = fill(NaN,nMax)
for n = 1:nMax
#local h, u, d, p, STree, Ce #only needed in REPL/script
h = m/n #time step size (in years)
u = exp(σ*sqrt(h))
d = exp(-σ*sqrt(h))
p = (exp(y*h) - d)/(u-d)
STree = BuildSTree(S,n,u,d)
Ce = EuOptionPrice(STree,K,y,h,p,false)
CeM[n] = Ce[1][] #pick out the call price at the starting node
end
In [14]:
plot( 1:nMax,[CeM fill(C,nMax)],
linecolor = [:red :blue],
linestyle = [:solid :dot],
linewidth = [1 2],
label = ["BOPM" "BS"],
title = "BOPM and BS call option price",
xlabel = "number of time steps (n)",
ylabel = "option price",
annotation = (170,2.7,text("h = m/n",8)) )
Out[14]:
In [ ]: