FMO6
FMO6 Lecture 3
Dr John Armstrong King’s College London August 3, 2016
FMO6
Matlab programming summary
We have learned how to use Matlab
What is the dierence between * and .*? Ifv = [8 -3 2 ]whatisv>0?
We have learned to write code as small functions
We have learned how to write unit tests
We have learned how to integrate using the rectangle rule
FMO6
Integration
Integration using rectangle rule
Algorithm (Rectangle Rule)
Let f : [a, b] −→ R be a function we wish to integrate. Choose a number of points N and dene:
h=b−a N
sm =a+(m−1)h 2
N
R = h f (sm) i=1
Then R is the rectangle rule estimate for b f (t ) dt . a
FMO6
Integration
Rectangle Rule
FMO6
Integration
Trapezium Rule
FMO6
Integration
Algorithm (Trapezium Rule)
Let f : [a, b] −→ R be a function we wish to integrate. Choose a number of panels n and dene:
h=
sm = a+mh
b−a n
T = h(f(s0)+2f(s1)+2f(s2)+…+2f(sn−1)+f(sn)) 2
Then T is the trapezium rule estimate for b f (t) dt. The number a
of integration points is N = n + 1.
FMO6
Integration
The next algorithm to consider is Simpson’s rule where one approximates the function as being piecewise quadratic.
Let f be a quadratic function then
hh f(t)dt =
(f(−h)+4f(0)+f(h))
−h
3
To see this:
Observe that both sides are linear in f .
Observe that for constant f , the result is true since 2 = 1 (1 + 4 + 1)
3
Observe that for odd functions, the formula will always hold.
In particular it holds for f (x ) = x . So by linearity it holds for all linear functions.
Observe that the result is true for f (x) = x2. So by linearity it holds for all quadratic functions.
More surprisingly the result is true if f is cubic. This follows from the fact that x3 is an odd function and by linearity.
FMO6
Integration
Algorithm (Simpson’s Rule)
Let f : [a, b] −→ R be a function we wish to integrate. Choose an even number of panels n and dene:
h=
sm = a+mh
b−a n
S = h(f(s0)+4f(s1)+2f(s2)+4f(s3)+2f(s4)+… 3
+ 2f (sn−2) + 4f (sn−1) + f (sn))
Then S is the Simpson’s rule estimate for b f (t) dt. The number a
of integration points is N = n + 1.
FMO6
Integration
Theorem
If f is four times dierentiable with |f (4)| < K for some K, then the error in the Simpson's rule estimate can be bounded above by
c, N4
where c is a constant depending upon a, b and K.
FMO6
Integration
By Taylor's theorem:
f(x) = f(a)+(x −a)f′(a)+ 1(x −a)2f(2)(a)+... 2
+ 1(x−a)3f(3)(a)+ 1(x−a)4f(4)(ξ) 3! 4!
for some ξ ∈ [a,x].
By construction Simpson's rule is exact for quadratic functions. Notice also that it is exact for cubic functions: this is because the integral of x3 over [−h,h] is zero as is the value given by Simpson's rule. By adding appropriate quadratic terms, any given cubic can be transformed to this special case.
The error over each interval of length 2h can thus be bounded by
h
error over interval ≤ c h4 = c h5
Thus the error over all intervals is bounded by c h4 = c . 3 N4
12 0
FMO6
Integration
Theorem
If f is twice dierentiable with |f (2)| < K for some K, then the error in the Trapezium rule estimate can be bounded above by
c N2
The same is true for the error in the Rectangle rule.
FMO6
Integration
Adaptive methods
Why distribute the points evenly? Maybe we'll get better convergence if we concentrate the points in areas where f is rapidly changing.
At stage n suppose we have grid Gn. We subdivide to get Gn+1. We can check to nd which splits were the most eective and use this to guide the next subdivision.
(Adaptive methods are not examinable)
FMO6
Integration
Gaussian integration
All our formulae have been of the form
b N
f≈ wif(xi)
a
Gauss asked: for polynomials of given degree d what are the best points and the best weights if we use d + 1 points. Gaussfound: (whena=0,b=1)choosethexi tobethe roots of the Legendre polynomials and the weights can be calculated in terms of the derivatives of the Legendre polynomials.
So for given d, Gauss gave a recipe to nd an integration rule which gives a perfect answer for polynomials of degree d with only d + 1 points.
(Gaussian quadrature is not examinable, sadly)
i=1
for some weights wi ∈ R and points xi ∈ [a, b].
FMO6
Integration
Practical methods
Matlab has a function integral which you can use in practice.
It uses an adaptive Gaussian quadrature strategy.
Similar functions can be found in other scientic libraries e.g. the GNU scientic library for C.
FMO6
Integration
The curse of dimensionality
By induction we can estimate d-dimensional integrals by
applying one of the rules above to an d − 1-dimensional
integral
f(x1,x2,...xd)dx1 dx2, ...dxd
=
f(x2,...xd)dx2, ...dxd dx1
Suppose we use N grid points for each one dimensional integration.
Let nd be the number of function evaluations required to compute a d-dimensional integral in this way.
n1 =N,nd =N∗nd−1. Thereforend =Nd.
e.g. for a 100-dimensional problem with 10 grid points we need 1 googol function evaluations (1 googol=10100).
FMO6
Integration
Throwing darts
If the area of the rectangle is 1m2
I estimate the area of the blob is 0.25m2
FMO6
Integration
Throwing darts
To nd the area of a given shape Σ
Bound it by a rectangle of known area A.
Throw n darts randomly (i.e. uniformly distribution) into the rectangle
Count the number of darts, nΣ, that land in the shape Σ. The area is approximately
nΣA n
FMO6
Integration
Monte Carlo integration
Algorithm (Monte Carlo Integration)
Let f : [a, b] −→ R. Estimate the integral of f by generating N random numbers xm which are uniformly distributed in the interval [a, b]. The Monte Carlo estimate for the integral is:
N
1 (b − a) f (xm)
N
i=1
Theorem
By the central limit theorem, under reasonable conditions on f , the
estimate converges at a rate of c 1 . N
√
FMO6
Integration
Multi-dimensional Monte Carlo integration
Algorithm (Monte Carlo Integration)
Let f : [0, 1]d −→ R. Estimate the integral of f by generating N random numbers xm which are uniformly distributed in the hypercube [0, 1]d . The Monte Carlo estimate for the integral is:
N
1 f(xm)
i=1
N
Theorem
By the central limit theorem, under reasonable conditions on f , the
estimate converges at a rate of O( 1 ). N
√
FMO6
Integration
Monte Carlo method
Experience suggests that dimension 3 or 4 is the approximate dimension where Monte Carlo integration begins to out-perform methods based on 1 dimensional integration rules.
Since you cannot generate true random numbers on a computer you have to make do with pseudo-random numbers. We will assume in this course that the pseudo-random number generator in MATLAB is up to the job.
In practice you should check how many random numbers you will need and make sure that you are using a generator that guarantees to provide high quality pseudo-randomness for that many numbers.
Monte Carlo doesn't completely get rid of curse of −1
dimensionality. The error is less than cN 2 but c can grow to be enormous for high dimensions.
FMO6
Integration
Indenite integrals
By performing a substitution one can transform an innite integral to an integral over an open interval.
Our integration rules (other than the rectangle rule) use the end points of the integral, so these limits had better exist after the substitution.
The choice of substitution makes a big dierence to the performance of the method. This doesn't just apply to innite integrals.
When integrating singular functions, it is a good practice to subtract the singularity from the integral. For example if the singularity can be approximated by the term 1/x1/2, subtract this o and integrate it separately analytically.
FMO6
Integration
The rectangle rule
function ret = integrateByRectangleRule( f, a, b, N )
h = (b-a)/N;
s = a + h*((1:N) - 0.5);
total = 0;
for x=s
total = total+f(x);
end
ret = h*total;
end
FMO6
Integration
The rectangle rule
function ret = integrateByTrapeziumRule( f, a, b, N )
n = N-1; % N = number of grid points, n=number of panels
h = (b-a)/n;
s = a + h*(1:n-1);
total = f(a)+f(b);
for x=s
total = total+2*f(x);
end
ret = 0.5*h*total;
end
FMO6
Integration
Simpson's rule
Exercise
FMO6
Integration
Indenite integration
function ret = integrateToInfinity( f, x, N, integrationRule)
% Performs a substitution to change
% the infinite integral to a finite integral.
if nargin<4
integrationRule=@integrateByRectangleRule;
end
function r = transformedFunction( s )
r = s^(-2) * f( x - 1 + 1/s );
if (isnan(r))
r = 0.0;
end end
ret=integrationRule( @transformedFunction, 0, 1, N );
end
FMO6
Integration
You can write a function that provides default values for arguments using nargin. This contains the number of parameters the user has actually passed in. So if it is lower than you expect, supply default values.
Use the function isnan to see if a calculation has evaluated to Not a number.
Our integrateToInfinity uses the rectangle rule by default and assumes that the limit of the transformed integral is 0.
FMO6
Integration
Calculating 1 sin(x ) dx 0
0 10
−2 10
−4 10
−6 10
−8 10
−10 10
−12 10
−14 10
−16 10
Errors in numerical integration
Rectangle rule
Trapezium rule Simpsons rule Monte Carlo
0123456 10 10 10 10 10 10 10
Number of points
Error
FMO6
Integration
Interpretation
The gradients of the log-log plot show the rate of the convergence.
Simpson's rule has gradient −4, Trapezium rule has gradient −2, Rectangle rule has gradient −2 as expected.
Monte Carlo is not inconsistent with the gradient −1 2
prediction. We could perform repeated runs to nd the average Monte Carlo error.
After a certain point the error of Simpson's rule increases. This is due to accumulating rounding errors. These grow with slope 1/2 consistent with the central limit theorem: we're adding independent errors.
Simpson's rule approaches machine precision accuracy.
FMO6
Integration
The code
function plotErrors()
% Plots the error of computing integral of sin from 0 to 1 points=1:18;
NValues = 2.^points + 1;
answer = -cos(2)+cos(0);
errorR = zeros( 1, length(NValues)); errorT = zeros( 1, length(NValues)); errorS = zeros( 1, length(NValues)); errorM = zeros( 1, length(NValues)); for i=1:length(NValues)
N = NValues(i);
fprintf('Running calculation with %d points\n', N);
errorR(i) = abs(integrateByRectangleRule(@sin,0,2,N) - answer); errorT(i) = abs(integrateByTrapeziumRule(@sin,0,2,N) - answer); errorS(i) = abs(integrateBySimpsonsRule(@sin,0,2,N) - answer); errorM(i) = abs(integrateByMonteCarlo(@sin,0,2,N) - answer);
end
figure();
loglog( NValues, errorR, NValues, errorT, NValues, errorS, NValues, errorM ); title('Errors in numerical integration');
xlabel('Number of points');
ylabel('Error');
legend('Rectangle rule', 'Trapezium rule', 'Simpsons rule', 'Monte Carlo' );
end
FMO6
Integration
Matlab functions for plotting
figure. Start a new drawing window. plot(x1,y1,x2,y2,...,xn,yn). Plot a line chart of the
vector yi against the vector xi . So this will contain n lines. loglog. Plot a log-log plot. Otherwise used just the same as
plot.
xlabel, ylabel, title. Self explanatory.
legend. Provide names for each series in the plot. For example:
legend('Series1','Series2','Series3');
FMO6
Integration
Question
8 How accurate do we need to be when pricing derivatives?
FMO6
Integration
How much accuracy?
There is a bid-ask spread on the stock price.
The model parameters used in a calculation is found by nding the best t to a the smile or the best t to historical data.
No need for machine precision typically.
Calculations of unlikely events or across entire portfolios with much hedging are sensitive to rounding errors.
FMO6
Integration
Risk neutral pricing
The risk neutral price of a derivative is its discounted expected value under a given risk neutral measure Q.
For the Black Scholes model you are given the price process in the physical measure P and can deduce the process in Q (by Girsanov's theorem)
For more complex models, one often simply chooses a Q measure process (by calibrating to the market) and ignores the P measure entirely.
FMO6
Integration
The Big Fight
Example
Wladimir Klitschko is due to ght Dr John Armstrong at the O2 Arena. A bookmaker is oering odds of 1000 to 1 on Armstrong winning. What can you say about the bookmaker's odds on Klitschko? What does this tell us about the probability of Armstrong winning?
FMO6
Integration
The Big Fight
Solution
The odds on Klitschko will be 1 to 1000 or less otherwise you can arbitrage against the bookmaker.
This tells us nothing about the probability of Armstrong winning, only about market prices for the bet. People may be willing to bet on Armstrong at 1000 to 1 even though the probability of him winning is probably somewhat less than this. After all you might win $1000 and it only costs $1.
The risk neutral probability of Armstrong winning is 1/1001 and the risk neutral probability of Armstrong losing is 1000/1001.
Q is determined by matching trades.
The bookmaker knows Q but has no opinion about P.
FMO6
Integration
Risk neutral pricing and integration
In risk neutral pricing the price is the discounted Q-expected payo of the derivative.
For a derivative with payo determined entirely by the stock price S at time T we have, by denition of expectation:
EQ (payo) =
Write q for the probability density function of S at time T.
Write P(S) for the payo given the nal stock price. price = e−rT EQ(payo)
−rT R
So we can price such a derivative by integration once we know the p.d.f. q.
R
payo(ST ) × Q-probability-density(S ) dS
= e
P(S)q(S)dS
FMO6
Integration
The Black Scholes Model, P measure
Measure P
PriceprocessSt: dSt =St(μdt+σdWt)
Ito's Lemma ⇓
Logpriceprocessst: dst =(μ−1σ2)dt+σdWt
2
Denition of s.d.e. and Wt ⇓
√ Distribution of s : N(s +(μ− 1σ2)T,σ T)
T0 Denition of ln N ⇓
2
√ Distribution of S : lnN(s +(μ− 1σ2)T,σ T)
T0 Substitute s = ln S
into of normal p.d.f. ⇓
2
(ln(S/S0)−(μ−1σ2)T)2
p.d.f.ofS : 1 exp − 2 T√2
Sσ 2πT 2Tσ MeanofST: exp(s0+μT)=S0exp(μT)
Integrate ⇓
FMO6
Integration
Log normal p.d.f.
To derive the p.d.f. of the log normal distribution, suppose that x ∼ N(α,β). So
t 1 (x−α)2
√ exp − −∞β 2π
P(x≤t) =
= √exp− dX
dx et 2
2β2
1 (ln X − α)
2β2
Here we've made the substitution x = ln(X ). We also have:
0 Xβ 2π
P(x ≤t)=P(ln(X)≤t)=P(X ≤et) So the p.d.f. of the log normal distribution is:
1 (lnX−α)2
√ exp − Xβ 2π
2β2
FMO6
Integration
The risk neutral measure
By Girsanov's theorem, if we can change the drift term in our price process by changing to an equivalent measure.
From our expression EP (ST ) = S0 exp(μT ) we see that the process with μ = r will be risk neutral. That is to say, if
μ = r, the expected stock price grows at the same rate as the zero coupon bond.
Therefore by taking μ = r we obtain the process in the risk neutral measure Q.
FMO6
Integration
The pricing kernel
Theorem
The Q-p.d.f. of ST is:
1
(ln(S/S0)−(r − 1σ2)T)2 2
2Tσ2
Given a European derivative that pays o P(S) if the stock price at
qT(S)= √ exp− Sσ 2πT
time T is S, the risk neutral price of this derivative is:
−rT ∞ 0
price = e
P(S)qT (S) dS.
Note that the derivative is assumed to be European and not path dependent. The function qT(S) is often called the pricing kernel.
FMO6
Integration
Pricing by numerical integration
Example
A dividend free stock follows the Black Scholes process with unknown drift and 10% volatility. The current stock price is $100. You should assume that the risk free rate is 5% APR. Use numerical integration to nd the risk neutral price of a 3 month European call option on the stock with strike $105?
FMO6
Integration
Solution - Units
In nance time is measured in years, so T = 0.25 years. Volatility is quoted as a percentage change over a year, so
−1 σ = 0.1years 2 .
The risk free rate, r, used in the Black Scholes Formula is a continuously compounded rate. So we have er = 1.05 if the annual percentage rate (APR) is 5%. Thus r = ln(1.05).
I'd expect most exam questions to say more simply σ = 0.1 and r = ln(1.05) ≈ 0.05, but you need to be more careful when using real market data.
FMO6
Integration
Solution - The callPayoff function
The payo of a European call option can be computed using the code below.
function p=callPayoff(K, S)
inMoney = S > K;
p = inMoney.*(S-K);
end
You could rewrite this using the maximum function if you wished.
FMO6
Integration
Solution – The pricingKernel function
We have computed the pricing kernel mathematically and found that it is:
1 qT(S)= √ exp−
(ln(S/S0)−(r − 1σ2)T)2 2
Sσ 2πT
we can translate this into MATLAB:
2Tσ2
function q=pricingKernel( S, T, S0, r, sigma )
coefficient = 1./(S * sigma * sqrt( 2 * pi * T));
numerator=-(log(S/S0)-(r-0.5*sigma^2)*T).^2;
denominator = 2*T*sigma^2 ;
q = coefficient .* exp( numerator/denominator);
end
FMO6
Integration
Solution – The priceCallByIntegration function
We can now use the theory of risk neutral pricing to compute the
payo.
−rT ∞ 0
price = e
P(S)qT (S) dS.
function p=priceCallByIntegration( K, T, S0, r, sigma )
function ret=integrand( S )
ret = callPayoff( K, S ) .*…
pricingKernel( S, T, S0, r, sigma );
end
p = exp(-r*T) …
* integrateToInfinity( @integrand, 0, 10001);
end
FMO6
Integration
How about some tests?
function testPriceCallByIntegration()
% A call option with strike 0 is the same thing
% as buying the stock. Check we get the correct answer.
% This is effectively a test that our pricing
% kernel is risk neutral
r = log(1.05);
S0 = 100;
sigma = 0.1;
T = 0.25;
actual = priceCallByIntegration(0,T, S0, r, sigma );
assertApproxEqual( actual, S0, 0.001 );
end
FMO6
Integration
Comparison with the Black Scholes Formula
function testBlackScholesCallPrice()
%Compares the black scholes price against
%that obtained by integration
r = log(1.05);
S0 = 100;
sigma = 0.1;
T = 0.25;
K = 110;
actual = priceCallByIntegration(K,T, S0, r, sigma );
expected = blackScholesCallPrice(…
K,T, S0, r, sigma );
assertApproxEqual( actual, expected, 0.001 );
end
FMO6
Integration
Now we’re condent, let’s answer the question
function answerQuestion()
%Let’s answer the question finally
T = 0.25;
r = log(1.05);
S0 = 100;
K = 105;
vol = 0.1;
answer = priceCallByIntegration(K,T,S0,r,vol);
fprintf(‘The answer is %f dollars\n’, answer );
fprintf(‘The answer to two d.p. is %.2f dollars\n’, answer );
end
This demonstrates the MATLAB code required to print things more attractively.
FMO6
Integration
What have we gained?
We can easily adapt the code to puts and digital options. We can price derivatives with arbitrary payos.
Given a pricing kernel qT(S) we can price derivatives. For example, you could take a pricing kernel with fat tails.
FMO6
Integration
A standard substitution
What substitution should we use to change our indenite integral to a denite integral? We have the following formula for the price:
−rT ∞
e payo(S)q(S)dS
0
where q(S) is the Q-p.d.f. Let Q(S) be the Q cumulative density function and make the substitution R = Q(S).
dR = q(S)dS
by denition of the probability density function as the derivative of
the c.d.f.. So the price is given by:
1 0
e −rT
payo(Q −1 (R ))dR
FMO6
Integration
Monte Carlo pricing
The price is given by:
exp−rT
So the monte carlo price is given by:
Generate random numbers R between 0 and 1.
Apply Q−1 to get a sequence of random numbers with the same distribution as the Q-distribution of stock prices.
Compute the average discounted price.
i.e. compute the expected value using a Q-measure simulation.
1 0
payo(Q −1 (R ))dR
FMO6
Integration
The Black Scholes case
Under the measure Q
S∼logN s0+ r−
σ2 √ T,σ T
2
So by denition of the log normal distribution
1 σ2
Q(S) = N √ ln(S/S0) − r − σT
Q−1(R) = S0 exp r −
We can compute N−1 using MATLAB’s norminv function.
T
2
σ2 √
T + σ TN−1(R)
2
FMO6
Integration
Exercises 1
8 Implement and test Simpson’s rule.
8 Use rand to generate uniformly distributed points in the unit square. Compute how many land in the unit circle. Hence estimate π.
8 Adapt the code to price a put option, a digital call option and a digital put option. Try using the OVME function on the Bloomberg terminals to test your answers.
8 Derive a formula for the delta of an option in terms of the derivatives of the pricing kernel. Use this to compute the delta of an option using integration.
FMO6
Integration
Exercises 2
8 Use the integral formulae for the price in terms of the payo and Q−1 to price a call option by Monte Carlo integration and using the rectangle rule.
8 Plot graphs of convergence of the two approaches.
8 Our estimates of the convergence of the integration methods assume a certain degree of dierentiability in the integrand. Split the integrals your are performing into two integrals at the singularity and see if that improves convergence.