ECS550NFB

Introduction to Numerical Methods using Matlab

Day 5
Lukas Laffers
[email protected]
Department of Mathematics, University of Matej Bel

June 12, 2015

Today

Derivatives pricing
I
I
I

Binomial lattices
Monte Carlo simulation
Solving a partial differential equation (e.g. finite differencing)

This lecture is based on chapters 2, 4, 7, 8 and 9 in

Brandimarte, Paolo. Numerical methods in finance and economics: a
MATLAB-based introduction. John Wiley & Sons, 2013.

Notation

I
I
I
I
I

C F P PV
rt -

coupon
face value
price
- present value
interest rate

PV =

T
X
t=0

FT
Ct
+
t
(1 + rt )
(1 + rT )T

cashFlows = [0 3 3 3 3 103];
interestRate = 0.03;
pvvar(cashFlows,interestRate)

Notation
I

- internal
rate of return irr(cashFlows), solution to
PT
FT
Ct
P = t=0 (1+)
t + (1+)T
or P
Ct /m
FT
P = Tt=0 (1+/m)
t + (1+/m)T if coupon is payed m times a year.
V (tn )tn
D - is duration, D = P V (t0 )t0 ++P
, where P V (ti ) is a
PV
present value of cashflow in time ti .
cfdur(cashFlows, interestRate),
DM - is modified duration, DM = D/(1 + /m). Why is it
interesting?
dP
d = DM P , so P DM P

C=
d2 P
d2

1 d2 P
P d2 ,

Why is it interesting?
= P C, so P DM P + P2C ()2

bndprice(Yield, CouponRate, Settle, Maturity)

Example: Portofolio with given duration and convexity

I
I

We are asked to create a portofolio with duration D and convexity

C from 3 bonds.
How will we set the weights?

P3

Pi=1
3

i=1
P
3

D i wi = D
C i wi = C

i=1 wi

(Note that this is an approximation)

=1

Mean-Variance portofolio optimization

I
I

We have n stocks with mean returns {r1 , r2 , rn } and covariance

matrix of their returns
How will we set our weights in the portofolio to get on average r
and minimizing risk (measured by the variance)?

minw wT w
s.t.
Pn

Pi=1
n

wi ri = r

i=1 wi

=1

wi 0
Makes sense if returns are normally distributied and the utility
function is quadratic.
[PRisk, PRoR, PWts]= frontcon(returns,covMatrix,100);

Derivatives

I
I

How do we find a right price?

How do we insure against potential losses?

What computational methods can we make use of

I
I
I

Solving a partial differential equation (e.g. finite differencing)

Monte Carlo simulation
Binomial lattices

Derivatives Pricing - simplest case

We construct a portofolio: 0 = S0 +
that replicate the option payoff
I
I

I
I
I

u = S0 u + ert = fu
d = S0 d + ert = fd

we get ,
by no-arbitrage principle f0 = 0
it is not the discounted expected value of the
payoff
if we use risk-neutral measure, it can be
interpreted as discounted expected value of
the payoff

Source: Brandimarte, p.105

Derivatives pricing - continuous time

Price of an Option, whose payoff depends only on the current price of

an underlying asset, or its price at maturity, follows Black-Scholes
differential equation (deterministic!).
f
f
1
2f
+ rS
+ 2 S 2 2 rf = 0
t
S 2
S
This PDE captures the dynamic, but we need to supply initial or
terminal conditions (f (S, T ) = max{S K, 0}).
In most cases, this equation must be solved numerically (but not for
vanilla European put/call options blsprice).

Derivatives pricing - binomial lattices

Source: Brandimarte, p.403

Derivatives pricing - binomial lattices

I
I
I

How do we calibrate the parameters of our binomial lattice? d, u, p

We can e.g. match the first two moments if St follows
dS = rSdt + SdW
Cox, Ross, Rubinstein (CRR): we have one degree of freedom, so
we canset u = 1/d and reduce the computational burden and get
rt d
u = e t , d = e t , p = e ud


Jarrow-Rudd: we set p =


d=e

2
r 2

t t

1
2

and calculate u = e

r 2

t+ t

and

Derivatives pricing - binomial lattices

Left: Stock, Right: European Call option. Source: Brandimarte, p.407

Binomial lattice - code

function [price, lattice] = LatticeEurCall(S0,K,r,T,sigma,N)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for i=0:N
lattice(i+1,N+1)=max(0 , S0*(u^i)*(d^(N-i)) - K);
end
for j=N-1:-1:0
for i=0:j
lattice(i+1,j+1) = exp(-r*deltaT) * ...
(p * lattice(i+2,j+2) + (1-p) * lattice(i+1,j+2));
end
end
price = lattice(1,1);

Binomial lattice - code improvement

Inefficiencies
I we repeat discounted probabilities exp(-r*deltaT)*p and
exp(-r*deltaT)*(1-p)
I we use a full matrix lattice - half of it is empty

Source: Brandimarte, p.412

Pricing American Options

I
I

We compare immediate payoff with continuation value

fi,j = max{K Si,j , ert (pfi+1,j+1 + (1 p)fi,j+1 )}

[AssetPrice, OptionValue] = binprice(Price, Strike, Rate,

Time, Increment, Volatility, Flag)
price = AmPutLattice(S0,K,r,T,sigma,N)

Derivatives Pricing - trinomial lattice

I
I

We will set pu , pm , pd to match the moments

of X, where Xt = log St

Choice of x? Rule of thumb x = t

Source: Brandimarte, p.422

Round-up: Should you use lattices?

I
I
I
I
I

Simple
Fast when problem dimensionality is small
Difficult for more complex path-dependent options
Discretization is important - e.g. poor choice of x may lead to
negative probabilities
Incorporating extensions usually requires careful implementation

Derivatives pricing - Monte Carlo Methods

I
I
I

simple, flexible
computationally intensive
when facing a complex problem, sometimes the only way to go

Monte Carlo Methods - Asset Paths

There are two sources of errors

I
I

sampling error
discretization error
Euler scheme

Brownian motion

Geometric Brownian motion

dSt = a(St , t)dt + b(St , t)dWt

St = a(St , t)t + b(St , t)) t

dSt = St dt + St dWt

St+t = (1 + t)St + St t

Discretizing Geometric Brownian Motion

But if we generate Geometric Brownian motion using

St+t = (1 + t)St + St t
St will be normal, rather than log-normal.
Solution? Using Itos

lemma on dSt = St dt + St dWt we get
d log St = 12 2 dt + dWt




Rt
which is integrated to St = S0 exp 12 2 t + 0 dW ( ) , so we
can generate sample paths using




1 2
St+t = St exp
t + t .
2
To speed things up, should we vectorize? In most cases, yes. The best
advice is to give it a try.

Example: Simple Monte Carlo for European Call

I
I

Create many sample paths for Stock price.

Discount and take average.

function [Price,CI] = BlsMC2(S0,K,r,T,sigma,NRepl)

nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T) ;
DiscPayoff = exp(-r*T)*max(0, S0*exp(nuT+siT*randn(NRepl,1))-K) ;
[Price,Var, CI] = normfit(DiscPayoff);
end

Monte Carlo Methods - Hedging strategies

I
I

Stop loss strategy - hold an asset if St > K, sell it if St < K

Delta hedging - replicate the option by holding stocks

Comparison of the costs of the strategy

HedgingScript.m
NSteps = 10
true price = 4.732837
cost of stop/loss (S)
cost of delta-hedging
NSteps = 1000
cost of stop/loss (S)
cost of delta-hedging

= 4.826756
= 4.736975
= 4.860783
= 4.732051

Exchange Option
At time T , we may exchange U for V and receive max{VT UT , 0}.
Suppose Ut and Vt follow bidimensional geometric Brownian motion
with drift r and the two Wiener processes have instantaneous
correlation
dU (t) = rU (t)dt + U U (t)dWU (t)
dV (t) = rV (t)dt + V V (t)dWV (t)
How to generate (correlated) sample paths?




1 p 0
1
T
,
=
, = LL , L =
1
1 2

1 = Z1
2 = Z1 +

p
1 2 Z2

[p,ci] = ExchangeMC(V0,U0,sigmaV,sigmaU,rho,T,r,NRepl)

Down-and-out Put option

Payoff is zero whenever the stock price hits the barrier level Sb .
I
I
I

Crude Monte Carlo

Conditional Monte Carlo
Importance sampling

Crude Monte Carlo - Down-and-out Put option

function [P,CI,NCrossed] = DOPutMC(S0,K,r,T,sigma,Sb,NSteps,NRepl)

Payoff = zeros(NRepl,1);
NCrossed = 0;
for i=1:NRepl
Path=AssetPaths1(S0,r,sigma,T,NSteps,1);
crossed = any(Path <= Sb);
if crossed == 0
Payoff(i) = max(0, K - Path(NSteps+1));
else
Payoff(i) = 0;
NCrossed = NCrossed + 1;
end
end
[P,aux,CI] = normfit( exp(-r*T) * Payoff);

Conditional Monte Carlo - Down-and-out Put option

I
I
I
I

Pdo = P Pdi
Discretization of the time interval of width t, T = M t, price
path is S = {S1 , S2 , . . . , SM }
Pdi = erT E[I(S)(K SM )+ ], where I(S) = 1 iff the price path hit
the barrier (Sj < Sb for some j)
If the barrier is crossed before maturity:

E[I(S)(K SM )+ |j , Sj ] = er(T t ) Bp (Sj , K, T t ), where j ,

where t = tj denotes the crossing time and Bp is the price of a

put option. payoff I(S)ert Bp (Sj , K, T t ), where j , where

t = t j .

The problem is that we hit the barrier too few times (249 out of
200000), in all other cases, the option price is zero.
[Pdo,CI,NCrossed] =
DOPutMCCond(S0,K,r,T,sigma,Sb,NSteps,NRepl)

Importance Sampling - Down-and-out Put option

We can change the sampling design, simulate paths that hit the barrier
more often and then correct for this change.
I

Eg

for the path sampling wegenerate

 independent normal draws Zj
2
with expected value = r 2 t and variance 2 t, then set
log Sj log Sj1 = Zj
instead, we will now generate from density with expected value
b and thus will cross the barrier more often
f (Z)I(S)(KSM )+
|j , Sj
g(Z)

The likelihood ratio

f (z1 ,...,zj )
g(z1 ,...,zj )

f (z1 ,...,zj )
g(z1 ,...,zj ) Ef [I(S)(K

SM )+ |j , Sj ]

is a correction term.

[Pdo,CI,NCrossed] =
DOPutMCCondIS(S0,K,r,T,sigma,Sb,NSteps,NRepl,bp)

Arithmetic Average Asian Option

The call option payoff is

(
max

N
1 X
S(ti ) K, 0
N

i=1

I
I

Crude Monte Carlo

Control Covariates
I
I

sum of asset prices as control covariate

geometric average option as a control covariate

Crude Monte Carlo - Arithmetic Average Asian Option

function [P,CI] = AsianMC(S0,K,r,T,sigma,NSamples,NRepl)

Payoff = zeros(NRepl,1);
for i=1:NRepl
Path=AssetPaths(S0,r,sigma,T,NSamples,1);
Payoff(i) = max(0, mean(Path(2:(NSamples+1))) - K);
end
[P,aux,CI] = normfit( exp(-r*T) * Payoff);

Control Variable - Arithmetic Average Asian Option

Sum of asset prices

PN

Y =

E[Y ]

i=0hS(ti )
i
PN
=E
i=0 S(ti )

r(N +1)t

= S(0) 1e1ert

[P,CI] = AsianMCCV(S0, K,r,T,sigma,NSamples,NRepl,NPilot)

Control Variable (2) - Arithmetic Average Asian Option

Geometric average option

I

nQ
o
Payoff is max ( S(ti ))1/N K, 0

Formula for the price of the geometric average option is available

[P,CI] = AsianMCGeoCV(S0,K,r,T,sigma,NSamples,NRepl,NPilot)

Estimating Greeks by Monte Carlo

Option sensitivities.
df (S0 )
dS0

f (S0 +S0 )f (S0 )

S0
E [C(S0 +S0 ,)][C(S0 ,)]
limS0 0
S0
[C(S0 S0 ,)]
it using limS0 0 E [C(S0 +S0 ,)]E
2S0

we can use

instead
or improve
and using common numbers
another approach is to calculate option price and in one
simulation - pathwise estimator

I
I

= limS0 0

Pathwise estimator - Greeks

I

discounted option payoff is a random variable

2
C = erT max{ST K, 0}, where ST = S0 e(r /2)T + T Z
C
dC dST
S0 = dST dS0

I dST
dS0
I dC
dST
I

=
=

ST
S0
erT I{ST

> K}

so the estimator for is erT SST0 I{ST > K}

function [Delta, CI] = BlsDeltaMCPath(S0,K,r,T,sigma,NRepl)

nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T);
VLogn = exp(nuT+siT*randn(NRepl,1));
SampleDelta = exp(-r*T) .* VLogn .* (S0*VLogn > K);
[Delta, dummy, CI] = normfit(SampleDelta);

Finite Difference Methods - Black Scholes formula

Assumptions
I
I
I
I
I

riskless rate r
stock price follows geometric Brownian motion with constant drift
and volatility
there are no dividends
no borrowing or lending constraints
no transaction costs

The the price of a derivative which payoff depends only on time and
stock price follows the Black-Scholes formula
f
1
2f
f
+ rS
+ 2 S 2 2 rf = 0
t
S 2
S

Finite Difference Methods

We will work with a discrete grid
I
I
I

S = 0, S, 2S, ..., M S
t = 0, t, 2t, ..., N t
fi,j = f (iS, jt)

There are different ways how to approximate partial derivatives

fi+1,j fi,j
f
S =
S
f f
I backward difference f = i,j i1,j
S
S
f
f
I central difference f = i+1,j i1,j
S
2S

forward difference

We also need boundary conditions

I
I

Call option: f (S, T ) = max{S K, 0}, S

Put option: f (S, T ) = max{K S, 0}, S

Boundary conditions

European vanilla put option

I
I
I

fi,N = max{K iS, 0}

f0,j = Ker(N j)t
fM,j = 0

European vanilla call option

I
I
I

fi,N = max{iS K, 0}
f0,j = 0,
fM,j = M S Ker(N j)t

Finite Difference Methods

Math preliminaries:
I
I
I

Explicit scheme - use forward approximation for the derivative

with respect to the time - may not be stable
Implicit scheme - use backward approximation for the derivative
with respect to the time - unconditionally stable
Crank-Nicholson method - combination of explicit and implicit
scheme - more precise approximation - unconditionally stable

Source: Brandimarte, chap. 5

Explicit scheme

fi,j fi,j1
t

fi+1,j fi1,j
2S
fi+1,j 2fi,j + fi1,j
1 22
i (S)2
= rfi,j
2
(S)2

+ riS
+

So the explicit scheme takes the following form

fi,j1 = afi1,j + bfi,j + cfi+1,j
price = EuPutExpl(SO,K,r,T,sigma,Smax.dS,dt)
Stability?

Stability - Explicit scheme

Explicit scheme is basically a trinomial lattice, hence we must choose
the discretization parameters S and t, so that the probabilities would
be positive.

Source: Brandimarte, p. 480

Implicit scheme

fi,j+1 fi,j
t

fi+1,j fi1,j
2S
fi+1,j 2fi,j + fi1,j
1 22
i (S)2
= rfi,j
2
(S)2

+ riS
+

So the implicit scheme takes the following form

afi1,j + bfi,j + cfi+1,j = fi,j+1
price = EuPutImpl(SO,K,r,T,sigma,Smax.dS,dt)
Stability?

Stability - Implicit scheme

Stability is not an issue at all. We can see that rewriting the system in
matrix form and inspecting the eigenvalues of the matrix, they are all
smaller than one in absolute value (Brandimarte p. 209).

Source: Brandimarte, p. 480

Crank-Nicolson scheme
We use derivatives in between the grid values
i
h 2
2
2
I f2 (Si , tj+1/2 ) = 1 f2 (Si , tj+1 ) + f2 (Si , tj ) + O(S 2 )
2 S
S
S
I f (Si , tj+1/2 )
t

f (Si ,tj+1 )f (Si ,tj )

t

+ O(t2 )

M1 fj1 = M2 fj
Pricing barrier put option using Crank-Nicholson scheme:
Boundary conditions:
I
I

f (Smax , t) = 0
f (Sb , t) = 0

price = DOPutCK(S0,K,r,T,sigma,Sb,Smax,dS,dt)

American Options

I
I
I

to account for free boundary is trivial in the Explicit scheme, this

is however subject to possible instabilities
we can use implicit scheme with iterative Gauss-Seidel method for
solving Ax = b
we can also use Crank-Nicholson scheme

price = AmPutCK(S0,K,r,T,sigma,Smax,dS,dt,omega,tol)

Literature

I
I
I
I
I
I
I

Brandimarte, Paolo. Numerical methods in finance and economics: a MATLAB-based introduction.

John Wiley & Sons, 2013.
chapter 2 - Financial Theory
chapter 4.5 Variance reduction techniques
chapter 7 - Option pricing by binomial and trinomial lattices
chapter 8 - Option pricing by Monte Carlo methods
chapter 9 - Option pricing by finite difference methods
Hull, John C. Options, futures, and other derivatives. Pearson Education India, 2006.

The End

Thank you for your attention!