FMO6
FMO6 Lecture 4
Dr John Armstrong King’s College London August 3, 2016
FMO6
Overview
What we’ve seen
We have learned how to program Matlab
We have learned the importance of testing
We have learned how to use numerical integration to perform risk neutral pricing
What’s to come Simulating prices for
Monte Carlo pricing (today) Risk management Simulation
Pricing using PDE methods Optimization
FMO6
Monte Carlo Pricing
Monte Carlo Pricing
FMO6
Monte Carlo Pricing
Simulating stock prices
We wish to simulate stock prices following the model dSt =St(μdt+σdWt)
at the time T.
Simulating stock prices means writing a function which
generates stock prices that have exactly the same probability
distribution as ST .
One way to simulate stock prices is rst to simulate the log of
the stock price sT instead, then exponentiate.
So we need to simulate sT following 1
dst=μ−σ2 dt+σdWt 2
By denition of Brownian motion this means generating sT
whicharenormallydistributedwithmeans +(μ−1σ2)T
and standard deviation σ
√
0
2
T.
FMO6
Monte Carlo Pricing
Algorithm (Simulating BlackScholes prices at time T)
Generate a normally distributed random number ε with mean 0 and standard deviation 1.
Set
Set
1√
s ̃ = s + μ − σ 2 T + σ T ε
2
T0
S ̃T =exp(s ̃T) The probability density of
S ̃T
is the same as that of
ST So S ̃T is a simulated stock price.
FMO6
Monte Carlo Pricing
Mnemonic
We have the process
Write:
1 dst=μ−σ2 dt+σdWt
2
dst →δsT =sT −s0
dt → δT = T − 0 √
dWt→ δTε 1√
s −s = μ− σ2 T+σ Tε
2
To get:
T0
FMO6
Monte Carlo Pricing
Note on the Mnemonic
The previous slide suggests approximating a continuous stochastic process with a discrete process with time intervals δT.
This idea does work and is called the Euler method. It approximately simulates a continuous stochastic process. Such methods improve as δT tends to zero.
Our algorithm for simulating stock prices exactly simulates Black Scholes stock prices no matter how large δT is.
The point is that we can integrate Brownian motion exactly, but not more general processes.
FMO6
Monte Carlo Pricing
Growth rate of noise
Why do we say that
proportional to proportionality.
√
√ dWt→ δTε?
The change in the stock price over time Nh is composed of N independent small changes over time periods of length h.
By the central limit theorem (assuming the small changes have nite variance), the change in the stock price over time Nh has
variance proportional to Nh.
Therefore, the cumulative eect of the noise grows at a rate
δT. Sigma is dened to be the constant of This is why sigma has units of years−1/2.
√
FMO6
Monte Carlo Pricing
Simulating at multiple time points
Simulating the stock price at times
0=t
knockedOut = 1;
end end
if (~knockedOut)
finalPrice = priceHistory( p, nSteps );
if (finalPrice>strike)
end end
end end
payoff(p)=finalPrice-strike;
FMO6
Monte Carlo Pricing
Comments
We don’t need to pass the number of paths and the number of steps to this function, we’ve used the size function to deduce that from the priceHistory
Since this code is completely repetitive across scenarios we know that we can vectorize it to improve eciency. We’ll do this shortly.
Lets test our code rst.
FMO6
Monte Carlo Pricing
Step 2 – Test the payo function
function testComputeKnockoutPayoff2Loops()
stockPrices = [100,101,102; 100,120,107; 100,103,108 ];
payoffs = computeKnockoutPayoff2Loops(105,110,stockPrices);
assertApproxEqual( payoffs(1), 0, 0.001);
assertApproxEqual( payoffs(2), 0, 0.001);
assertApproxEqual( payoffs(3), 3, 0.001);
end
FMO6
Monte Carlo Pricing
Step 3 – Write the pricing function
function [price, errorEstimate]=priceKnockoutByMonteCarlo(…
strike, barrier, T,…
S0, r, sigma, …
nPaths, nSteps )
% Generate paths in risk neutral measure (mu=r)
priceHistory = generateBSPaths(T,S0,r,sigma,nPaths,nSteps);
payoffs = computeKnockoutPayoff(strike,barrier,priceHistory);
discountedPayoff = exp(-r*T)*payoffs;
price = mean( discountedPayoff );
errorEstimate = std( discountedPayoff )/sqrt(nPaths);
end
FMO6
Monte Carlo Pricing
Step 4 – Test the pricing function
Exercise
FMO6
Monte Carlo Pricing
Remarks
The main pricing function is remarkably simple
The main pricing function does little more than explain in English what the Monte Carlo pricing algorithm actually is.
The real work is done in generateBSPaths and computeKnockoutPayoff.
In the exam I may ask you to write pseudo-code for certain algorithms. I strongly recommend splitting your pseudo code into functions that explain the algortihm. Given how simple well written code is, you might well prefer to write actual code.
We’ll now see how to vectorize the computeKnockoutPayoff function.
FMO6
Monte Carlo Pricing
Removing the loop across paths
function [ payoff ] = computeKnockoutPayoff1Loop( …
strike, barrier, priceHistory )
nPaths = size( priceHistory, 1 );
nSteps = size( priceHistory, 2 );
knockedOut = zeros( nPaths, 1);
for t=1:nSteps
knockedOutThisTime = (priceHistory(:,t) > barrier);
knockedOut = knockedOut | knockedOutThisTime;
end
finalPrice = priceHistory( :, nSteps );
inMoney = finalPrice>strike;
payoff=(~knockedOut).*inMoney.*(finalPrice-strike);
end
FMO6
Monte Carlo Pricing
Understanding the vectorized code
We have used the vectorization trick of using the fact that knockedOut and inMoney take the values 1 and 0 since 1 represents true and 0 represents false.
This means
(~knockedOut).*inMoney.*(finalPrice-strike)
Will only be non zero when we haven’t knocked out and aren’t in the money.
We have used | which means element-by-element ‘OR’ just as .* means element-by-element multiplication.
FMO6
Monte Carlo Pricing
Vectorizing again
We are currently looping to nd out if the price was ever above the barrier.
Suppose we compute priceHistory > barrier. This will consist of 1’s when the price is above the barrier and 0’s when the price is below the barrier.
This means that the maximum of priceHistory>barrier in each row will be 1 if the barrier knocked out and 0 otherwise.
You can use max(A,[],2) to compute a vector of the maximum across the columns. max(A,[],1) computes the maximum across the rows.
FMO6
Monte Carlo Pricing
No loops
function [ payoff ] = computeKnockoutPayoff( …
strike, barrier, priceHistory )
knockedOut = max( priceHistory>barrier, [], 2);
notKnockedOut = 1-knockedOut;
finalPrice = priceHistory(:,end);
inMoney = finalPrice>strike;
payoff = inMoney .* notKnockedOut .* (finalPrice-strike);
end
FMO6
Monte Carlo Pricing
Which version to learn?
I’ve shown how to price the code using three short functions computeKnockoutPayoff, priceKnockoutCallOption and generateBSPaths.
We have three equivalent versions of generateBSPaths and computeKnockoutPayoff
I would focus on understanding the nal versions of the code. I introduced the other versions to help you understand the nal version.
FMO6
Monte Carlo Pricing
Should you vectorize your own code?
Good things about vectorization
Vectorization makes the code more readable once you know the standard tricks.
Vectorization may make the code faster Bad things about vectorization
Vectorization makes the code less readable if you don’t know the standard tricks.
Vectorization may make the code take longer for you to write (until you’ve mastered it)
I recommend writing what you nd comes most naturally and optimizing only if there is a problem (or the question tells you to)
FMO6
Monte Carlo Pricing
Unit testing and random numbers
Algorithms that use random numbers will sometimes give poor answers by chance.
If you’re not careful this means that your tests will sometimes fail by chance.
To make your tests reliable you should seed the random number generator.
Pseudo random number generators have a state which determines what they will do next.
Each time you generate a number they change state. This ensures they will generate a dierent number the next time. Manually xing the state is called seeding the random number generator.
FMO6
Monte Carlo Pricing
Seeding the random number generator
In MATLAB rng(‘default’) sets the random number generator back to its default state.
Therefore you should start unit tests of functions that use the random generator with a call of rng(‘default’).
FMO6
Monte Carlo Pricing
Computing Greeks by Monte Carlo
A trader will not thank you if you can compute the price but not the Greeks. Why is this?
FMO6
Monte Carlo Pricing
Computing Greeks by Monte Carlo
A trader is likely to be using some variant of the delta hedging strategy so they will want to make sure their portfolio is delta neutral
If they can’t ensure their portfolio is delta neutral, they won’t trade because they can’t hedge the risk of their position.
In general, a price is useless without a trading strategy to achieve the price.
FMO6
Monte Carlo Pricing
Numerical dierentiation
Let f be a smooth function.
We can approximate the derivative at x as
f ′(x) = f (x + h) − f (x) h
for small h.
A better approximation is
f ′(x) = f (x + h) − f (x − h) 2h
In the latter case, the rst order error terms cancel and the error is bounded by
sup
|f(3)(c)| 2 h
c ∈[x −h,x +h]
6
FMO6
Monte Carlo Pricing
What h should you choose?
We want to choose h that is small enough to give a reasonably accurate value value.
We want to choose h so it isn’t so small that rounding errors dominate the calculation.
√
Choosing h >
computer will ensure that h isn’t too small. For our purposes ε = 2.2 × 10−16
For ultimate accuracy one worries about whether h + x can be represented accurately on the computer and tweaks the value accordingly.
Since we will apply this to a Monte Carlo method, we needn’t worry about choosing the optimal h.
εx where ε denotes the accuracy of the
FMO6
Monte Carlo Pricing
How not to compute Monte Carlo delta
Do not:
Compute the Monte Carlo price with an initial stock price of S0
Compute the Monte Carlo price with an initial stock price of S0 + h
Take the dierence and divide by h
This is because the Monte Carlo prices are randomly generated. The random error will overwhelm the systematic dierence we are trying to measure.
FMO6
Monte Carlo Pricing
How to compute Monte Carlo delta
Algorithm (Monte Carlo Delta)
Compute the Monte Carlo price with an initial stock price of S0 − h
Compute the Monte Carlo price with an initial stock price of S0 + h using exactly the same random numbers in the simulation
Take the dierence and divide by 2h
FMO6
Monte Carlo Pricing
Why this works
Monte Carlo pricing is really just numerical integration
Under reasonable conditions, the partial derivative of an integral is the integral of the partial derivative.
Our Monte Carlo Delta algorithm is equivalent to computing the derivative of the pricing kernel numerically and then integrating by Monte Carlo integration.
FMO6
Monte Carlo Pricing
MATLAB implementation
function [ delta ] = computeDeltaByMonteCarlo( …
strike, barrier, T,…
S0, r, sigma, …
nPaths, nSteps )
h = 10^(-6)*S0; % Won’t cause rounding problems
% but a minute change financially
rng(‘default’);
p1 = priceKnockoutByMonteCarlo(strike,barrier,T,…
S0-h, r, sigma, …
nPaths, nSteps );
rng(‘default’);
p2 = priceKnockoutByMonteCarlo(strike,barrier,T,…
S0+h, r, sigma, …
nPaths, nSteps );
delta = (p2-p1)/(2*h);
end
FMO6
Monte Carlo Pricing
Explanation
We have seeded the random number generator to ensure that the same random numbers are used at S0 − h and at S0 + h.
It would be better in practice to avoid generating the same random numbers twice as this inevitably be faster and would get rid of the bias of always using the same random numbers in calculations.
FMO6
Monte Carlo Pricing
Second derivatives
f ′(x − h) ≈ f (x) − f (x − h) h
f ′(x) ≈ f (x + h) − f (x) h
f′′(x)≈f′′(x−h)≈ f′(x)−f′(x−h) h
≈ f(x +h)−2f(x)+f(x −h) h2
Somewhat more rigorously, one can remark that this formula for f ′′ is accurate for quadratics and so, using Taylor’s theorem one can give bounds on the error of the approximation.
FMO6
Monte Carlo Pricing
Antithetic sampling
Let P1 and P2 be random variables with E(P1) = E(P2). We wish to calculate this expectation.
Dene
Then
Z = P1 + P2 . 2
E(Z) = E(P1) = E(P2)
Var(Z)= Var(P1)+Var(P2)+2Cov(P1,P2)
4
So, all other things being equal, if P1 and P2 are negatively correlated, the estimate E(Z) will be improved over the case when P1 and P2 are independent.
FMO6
Monte Carlo Pricing
Monte Carlo with antithetic sampling
Suppose that in our Monte Carlo pricing algorithm, we generate a random normally distributed vector ε and use this to compute a payo P1 = P(S(ε)). We can use the vector −ε to compute another payo P2 = P(S(−ε)). Here S is the stock price given ε, P then computes the payo.
Cov(P1,P2) will be negative for many payo functions (if ε leads to a good payo, −ε will often be bad)
So the estimate E(Z) may be a better estimate than we would have obtained by taking independent samples P1.
So by computing the Monte Carlo price using a sample of random numbers εi and the antithetic sample −εi we will get a better estimate than we would obtain by simply doubling the sample size.
FMO6
Monte Carlo Pricing
Exercises (p1)
8 What tests could you write for the generateBSPaths function? Write them and check that they pass. Make sure that the tests always pass whatever numbers are generated.
8 Use the generateBSPaths function to create a plot of the percentiles of the stock price over time as shown in the slides.
8 Sketch a graph of the how the price of a knockout option varies as the barrier changes. What is the minimum value and when is it obtained? What is the maximum value and when is it obtained? Explain your answer. Check that your sketch matches the results of the Monte Carlo simulation.
FMO6
Monte Carlo Pricing
Exercises cont…
8 Price a knock in option and an Asian option. How can you use function passing to improve your code?
8 For our Monte Carlo pricing algorithm we need to generate normally distributed random numbers. Suppose that we do this by generating a uniformly distributed random number and then applying N−1. Show that the Monte Carlo pricing algorithm then becomes equivalent to one of the Monte Carlo integration techniques from the last lecture.
8 Refactor the code so that the random numbers used to compute the delta are not generated twice.
FMO6
Monte Carlo Pricing
Exercises cont…
8 We have stated an approximation formula for the second derivative in terms of the value of f at the points {x − h, x, x + h}. Assuming f is smooth, compute a bound on the error of this formula.
8 Use the approximation formula for the second derivative to write a function that computes the gamma of a European option using a Monte Carlo method. Check your answer.
8 Numerically compute the derivative of sin(x) in two dierent ways and plot a log-log plot of the error in the estimates.
8 Implement Monte Carlo pricing with antithetic sampling for a knockout option.