程序代写代做代考 matlab scheme algorithm FMO6

FMO6
FMO6 Lecture 6
Dr John Armstrong King’s College London August 3, 2016

FMO6
Simulating more interesting stochastic processes
Simulating more interesting stochastic processes

FMO6
Simulating more interesting stochastic processes
Summary so far
Simulating the Black􏺋Scholes model has already given some interesting results
Simulating in the -measure allows us to price derivatives
Simulating in the -measure allows us to test trading strategies

FMO6
Generating correlated random variables
Generating correlated random variables

FMO6
Generating correlated random variables
Correlated normally distributed random variables

FMO6
Generating correlated random variables
Multivariate normal distribution
Independent, normally distributed random variables with equal variance:
22 p(x,y) = ce−1x2e−1y2
= ce−1(x2+y2) 2
Where c is a normalization constant chosen so that 􏶩􏶩
p=1
Lines where x2 + y2 = constant.

FMO6
Generating correlated random variables
Di􏺈erent variances
−1(x )2 −1(2y)2 p(x,y)=ce 2 2 e 2
= ce−1(1×2+4y2) 24
Lines where
2
+ 4y = constant.
x2 2

FMO6
Generating correlated random variables
Correlated
p(x,y) = ce−1(3×2−2xy+3y2) 2
Lines where 3×2 − 2xy + 3y2 = constant.

FMO6
Generating correlated random variables
Conclusion
The contour lines of the p.d.f. of a multivariate normal distribution are ellipses.
Proof.
The p.d.f. is
2 p(x,y) = ce−1(q(x,y)
for some quadratic in x and y and normalization constant c. The curves p(x,y) = c1 are the contour lines of the p.d.f. These are the lines
􏰒c1􏰓 q(x,y)=−2log c =c2
This is true in higher dimensions too. Note that in fact, q must be positive semi-de􏺉nite by the de􏺉nition of a multivariate normal distribution, so these really are ellipsoids.

FMO6
Generating correlated random variables
Ellipse to distribution
Corollary
Given an ellipse E there is a unique 2 − d multivariate normal distribution with contour line E containing half the probability mass.

FMO6
Generating correlated random variables
Dictionary

FMO6
Generating correlated random variables
Dictionary quadratics and ellipses
We 􏺆recall􏺊 that quadratics correspond to conic sections. positive de􏺉nite quadratics in x and y describe an ellipse
De􏺉nition
(Rather lazy) A quadratic in x and y is positive de􏺉nite if it de􏺉nes an ellipse. We’ll give a better, algebraic de􏺉nition later.

FMO6
Generating correlated random variables
Summary
We have established a dictionary between 3 di􏺈erent mathematical objects:
Positive de􏺉nite quadratics q(x,y)
Ellipses q(x,y) = c. (The constant c in our dictionary isn’t interesting – it just corresponds to the choice of 50% of the probability mass inside the ellipse.)
2-dimensional multivariate normal distributions
ce−1q(x,y) 2
The normalization constant isn’t interesting.

FMO6
Generating correlated random variables
Matrix notation
The general quadratic is:
ax2 +2bxy +cy2 +dx +ey +f
Which can be written in matrix form as:
􏰉 􏰈􏰔ab􏰕􏰔x􏰕􏰉 􏰈􏰔x􏰕 xybcy+dey+f
Or equivalently:
Where
⃗x T Ω − 1 ⃗x + ⃗v T ⃗x + f
􏰔a b􏰕 􏰔d􏰕 􏰔x􏰕 Ω−1= b c , ⃗v= e , ⃗x= y

FMO6
Generating correlated random variables
N-dimensional dictionary
We have the following dictionary in higher dimensions: Positive de􏺉nite quadratic form:
⃗x Ω − 1 ⃗x + ⃗v T ⃗x with Ω a symmetric matrix
Ellipsoids for some (uninteresting) constant c1 depending only on the dimension.
Multivariate normal distributions: for some (uninteresting) constant c2 depending only on the dimension.
Note that Ω is the covariance matrix.

FMO6
Generating correlated random variables
Conclusion
We can use the geometry of ellipses to understand multivariate normal distributions
We can use the algebra of quadratic forms to provide formal proofs of results about normal distributions.

FMO6
Generating correlated random variables
Idea
Let X, Y be variables with a 2-d multivariate normal distribution with a given covariance matrix. We want to simulate these variables.
Let E be the corresponding ellipse.
Find coordinates x ̃ and y ̃ such that E is the unit circle:
􏰔 x 􏰕 􏰔 x ̃ 􏰕
y =A y ̃ +B
Simulate variables X ̃ and Y ̃ using randn
Then
􏰔 X ̃ 􏰕
A Y ̃ +B
will simulate the desired distribution

FMO6
Generating correlated random variables
Transformation
−→
( x , y ) ← → ( x ̃ , y ̃ )

FMO6
Generating correlated random variables
Option 1

FMO6
Generating correlated random variables
Option 2

FMO6
Generating correlated random variables
Algebraic interpretation
Option 1: Find the principal axes and rotate. This is equivalent to diagonalizing the symmetric matrix Ω
Option 2: Sheer and then rescale. This is equivalent to 􏺉nding the Cholesky decomposition
It’s much less algebraically tedious.

FMO6
Generating correlated random variables
Another interpretation
To simulate correlated variables X and Y First simulate X
Now compute Y0 = αX for some appropriate α determined by covariance
Now simulate D a normally distributed variable representing a random deviation from perfect correlation. Take
Y = Y0 + D = αX + D
This is essentially equivalent to Cholesky decomposition.

FMO6
Generating correlated random variables
Algebraic de􏺉nitions
De􏺉nition
A symmetric matrix A is positive de􏺉nite/semi-de􏺉nite if x T Ax > 0 (or geq0) for all vectors x.

FMO6
Generating correlated random variables
Algebraic version
We are given a quadratic form XTΩ−1X with Ω a symmetric and positive de􏺉nite n × n matrix and X an n vector.
We wish to 􏺉nd a matrix A such that if X = A X ̃
then the quadratic form is a circle in X ̃ coordinates. XTΩ−1X =X ̃ATΩ−1AX
So we require:
AT Ω−1A = 1 equivalently Ω−1 = (A−1)T A−1 equivalently AAT = Ω

FMO6
Generating correlated random variables
Pseudo square root
De􏺉nition
If Ω is a positive de􏺉nite symmetric matrix, then a matrix satisfying AAT =Ω
is called a pseudo square root of Ω.
Lemma
Given a pseudo square root A of Ω then if X ̃ is a vector of independent normally distributed random variables with mean 0 and standard deviation 1 then AX ̃ is a multivariate normal variable with mean 0 and covariance Ω.

FMO6
Generating correlated random variables
Square root by diagonalizing
Let Ω be positive semi-de􏺉nite and let e1, e2, . . . en be an orthonormal basis of eigenvectors with eigenvalues λ1, . . . λn then the matrix
A=􏰉√λe √λe … √λe 􏰈 1122nn
is a pseudo square root of Ω. Note that each entry shown in A is in fact a complete column vector.
PROOF:
De􏺉ne
 1 eT 
λ1 1
 1 eT  √
 λ2 2  B = ( A ) =  . . . 
1 eT √λn n
Note that we have been able to simply write down the inverse of A since the basis e1, e2 …en is orthonormal so eiTej = δi,j.
√
−1 T

FMO6
Generating correlated random variables
Proof continued
So we can write down B
􏰒111􏰓 B= √ e1 √ e2 …√ en
λ1 λ2 λn By de􏺉nition of an eigenvector Ωei = λi ei .
So
 1 eT 
λ1 1
1T √e􏰉√√√􏰈
√
T
 λ2 2 
BΩB=  λ1e1 λ2e2…λnen
… 1 eT
√λn n

FMO6
Generating correlated random variables
Cholesky Decomposition
Theorem
If Ω is a symmetric positive semi-de􏺉nite matrix then there exists a unique lower triangular matrix L with Ω = LLT and positive diagonal.
De􏺉nition
L is called the Cholesky decomposition of Ω.
Because L is lower triangular, we can write the equation LLT = Ω out in detail as:
a 0 0…0a a a …a 11 11 21 31 n1
 a21 a22 0 … 0  0 a22 a32 … an2  a a a …00 0a …a
31 32 33  33 n3=Ω  . . . .  . . . .  … .… .
an1 an2 an3 … ann 0 0 0 … ann

FMO6
Generating correlated random variables
Proof continued
n(n−1)
Take the lower triangular part of both sides. This gives equations in the same number of unknowns. The 􏺉rst row gives:
a2 =Ω 11 11
We can now solve for a unique positive a11. The next row gives: a21a11 = Ω21
a2 +a2 =Ω 21 22 22
We solve the 􏺉rst for a12. Now we can read o􏺈 the unique positive a22. The third row gives:
a31a11 = Ω31 a31a21 + a32a22 = Ω32
a2 +a2 +a2 =Ω 31 32 33 33
2

FMO6
Generating correlated random variables
Proof continued
We solve for a31 then a32 then a33.
Proceeding in this way gives an algorithm for computing the Cholesky decomposition.
To compute all the aij will take O(i2) computations for each i (this is the number of coe􏺍cients in the equations we write down). So the algorithm will take O(n3) steps.
Note that a complete proof requires additionally showing that that Ω being positive de􏺉nite means the quadratics we solve have real solutions, we’ll skip this detail.

FMO6
Generating correlated random variables
Summary
If we can 􏺉nd a pseudo square root A of a covariance matrix Ω we can generate normally distributed random numbers with coariance matrix Ω by [􏺉ll in the blanks].
One way to 􏺉nd a pseudo square root is be diagonalizing the matrix Ω (tedious)
Another way is through Cholesky decomposition (less tedious)

FMO6
Generating correlated random variables
Exercises
8 Use matlab’s chol function to 􏺉nd the Cholseky decomposition
of
􏰔3 1􏰕 12
8 Use matlab to plot a scatter plot of 10000 points (X , Y ) where X and Y are normally distributed with covariance matrix
􏰔3 1􏰕 12

FMO6
Generating correlated random variables
Exercises continued
8 X is normally distributed with mean 5 and standard deviation 3 and Y is normally distributed with mean 7 and standard deviation 1 and if X and Y have correlation ρ = 0.5. Generate a sample of points (X , Y ) matching these properties. How have you tested your answer?
8 (Harder) Geometrically Cholesky decomposition proceeds by performing a sheer in the direction of the axes x1, x2, x3, . . . in turn. If you change the order of the axes, you will get a di􏺈erent pseudo square root. Use this idea to 􏺉nd another pseudo square root of the matrix:
511 161 114

FMO6
Generating correlated random variables
Exercises continued
8 Write a function randnMultivariate(omega,n) which generates n samples from a multivariate normal distribution with covariance matrix omega.
8 Compute the Cholesky decomposition of
by hand.
􏰔3 1􏰕 12

FMO6
Solving SDEs numerically
Solving SDEs numerically

FMO6
Solving SDEs numerically
Reference
Kloeden and Platen 􏺆Numerical Solution of Stochastic Di􏺈erential Equations”.
We won’t give proofs.

FMO6
Solving SDEs numerically
Problem setting
1 dimensional stochastic di􏺈erential equation dXt =a(X,t)dt+b(X,t)dWt
n dimensional stochastic di􏺈erential equation dXt =a(X,t)dt+b(X,t)dWt
where
is an n×d matrix
Wr
is a d-dimensional Brownian motion
In either case we also have an initial condition X0.
Xt ∈Rn a(X,t) ∈ Rn b(X,t)

FMO6
Solving SDEs numerically
Notation
We will choose a time step δt and will 􏺉nd approximate solutions by discretization.
We will write di􏺈erence equations for our approximations. Time point i corresponds to the time iδt.
If Wt is a Brownian motion that we have been given then we de􏺉ne
δWi = Wiδt − W(i−1)δt so δWi is a random variable.

FMO6
Solving SDEs numerically
Euler scheme
The Euler scheme for solving the SDE is to de􏺉ne
X ̃i =X ̃i−1+a(Xi−1,t)δt+b(Xi−1,t)δWi
so each X ̃i is a random variable determined by the δWj with j ≤ i. We claim that (in a sense to be explained later) the X ̃i are a good approximation to Xiδt.

FMO6
Solving SDEs numerically
Application
To simulate the values Xi , generate independently normally distributed εi with mean 0 and standard deviation 1. De􏺉ne:
̃ ̃√ Xi = Xi−1 + a(Xi−1, t)δt + b(Xi−1, t)(
√
δt)εi
δt!
Solving the SDE when we are given values of Wt over time.
NOTE THE
Note the slight distinction between:
Simulating the stochastic process, where we run many simulations and generate our own values of εi .

FMO6
Solving SDEs numerically
Theorem
Suppose that:
|a(x, t) − a(y, t)| + |b(x, t) − b(y, t)| < K1|x − y| and and |a(x, t)| + |b(x, t)| < K2(1 + |x|) 2 |a(x,s)−a(x,t)|+|b(x,s)−b(x,t)| ξ2 to keep volatility positive.

FMO6
Solving SDEs numerically
r = 0 so we require that S is a martingale for this to be a valid Q measure model. This is why there is no drift term.
In Black􏺋Scholes model there is a unique compatible Q model for a given P. This isn’t true in general, so one usually takes Q as given.

FMO6
Solving SDEs numerically
The volatility smile
If we simulate Q measure stock prices in the Heston model and use this to compute risk neutral prices for options, will this 􏺆explain􏺊 the volatility smile?
(Question: what is implied volatility? What is the volatility smile?)

FMO6
Solving SDEs numerically
Monte carlo pricing code
Our pricing code looks like this:
function ret = priceByMonteCarlo (…)
paths =
payoffs = computeOptionPayoffs (…); ret = mean(payoffs)*exp(-r*T);
end
generatePricePaths (…);
All we need to do is to change this to use the Heston model to generate price paths.

FMO6
Solving SDEs numerically
function [ prices, variances ] = generatePricePathsHeston( …
S0, v0, …
kappa, theta, xi, rho, …
T, nPaths, nSteps)
%GENERATEPRICEPATHSHESTON Generate price paths according to the
% Heston model
prices = zeros( nPaths, nSteps );
variances = zeros( nPaths, nSteps );
currS = S0;
currv = v0;
dt = T/nSteps;
for i=1:nSteps
epsilon = randnMultivariate( [1 rho; rho 1], nPaths );
dW1 = epsilon(1,:)*sqrt(dt);
dW2 = epsilon(2,:)*sqrt(dt);
currS = currS + sqrt( currv).* currS .* dW1′;
currv = currv + kappa*(theta – currv)*dt + xi*sqrt( currv).* dW2′;
currv = abs( currv ); % Forcibly prevent negative variances
prices( :, i) = currS;
variances( :, i) = currv;
end

FMO6
Solving SDEs numerically
I ran the simulation with the following parameters
nPaths=100000 nSteps=50 T=1
S0 = 1
κ=2
θ = 0.04 v0 = 0.04 ρ=0
I then used three di􏺈erent valued of ξ.
I plotted the Black􏺋Scholes implied volatility for a number of strikes

FMO6
Solving SDEs numerically
0.21 0.208 0.206 0.204 0.202
0.2 0.198 0.196 0.194
Heston model volatility smile
xi=0
xi=0.2
xi=0.4
95% confidence, xi=0 95% confidence, xi=0
0.192
80 85 90 95 100 105 110 115 120

FMO6
Solving SDEs numerically
When ξ = 0 this becomes equivalent to the Black Scholes
model. So in theory the implied volatility equal to 0.2.
As well as computing monte carlo prices, bounds for the computation when ξ = 0. prediction 􏺉ts within the error bounds as
Other valued of ξ do give rise to a smile.
should be a constant
I’ve plotted error The Black Scholes
one would hope.

FMO6
Solving SDEs numerically
Other applications
Repeat the delta hedging simulation we did last week, but generate stock prices using di􏺈erent -measure models. We can use this to see how well delta hedging performs with e.g. fat tailed stock prices.
Repeat the VaR Monte Carlo simulations we’ll carry out next week with more interesting models
etc. etc.

FMO6
Solving SDEs numerically
The Milstein Scheme
The Euler scheme isn’t the end of the story. The Milstein scheme for:
dXt =a(Xt,t)dt+b(Xt,t)dWt
is to take
X ̃ i = X ̃ i − 1 + a δ t + b δ W t + 1 b ∂ b 􏰉 ( δ W t ) 2 − δ t 􏰈
2 ∂x This is Euler scheme plus one more term
Under certain bounds on the coe􏺍cients, converges in expectation at rate O(δt)
n-d versions exist but are more complex.

FMO6
Solving SDEs numerically
Plot of errors of Euler and Milstein
−1 10
Accuracy of schemes
Euler error
Milstein Error
−2 10
−3 10
01234 10 10 10 10 10

FMO6
Solving SDEs numerically
Exercises
8 Simulate the process
dSt = St(μdt + σdWt)
using the Euler scheme and 􏺉nd the exact solution too. Use this to generate a log-log plot of errors for the Euler scheme.
8 Simulate the Vasicek interest rate model drt =a(b−rt)dt+σdWt
Generate plots of interest rate paths with varying parameters so you get a feel for this kind of model

FMO6
Solving SDEs numerically
8 Modify the delta hedging code from last week so that one still delta hedges as though one believed the Black􏺋Scholes model was true, but in fact the interest rates are stochastic and follow the Vasicek model. How does the delta hedging strategy perform?