Computer practical∗ Svetlana Borovkova Artem Tsvetkov
FEWEB / VU FEWEB / VU September 2, 2016
1 Matlab refreshment
Please refresh your Matalb knowledge before the practice. We will move quickly through this in the class. Please look through the accompanying m-file. There is more on Matlab programming in it.
HELP is very useful in Matlab. Access it through Search Documentation in the upper right corner and learn how to use it.
1.1 Comments
Matlab ignores everything on the line after the percentage symbol %. We will use it for commenting the code.
1.2 Run the code
You can run Matlab code in different ways.
• Press F5, if want to run the code from the beginning till the end,
• Or in Editor tab click Run.
• To run a section from double comment, %%, to the next double comment, click anywhere inside and press Ctrl-Enter,
• Or click Run Section in the Editor menu.
1.3 Vectors, matrices, and operations with them
Matlab is a natural tool to operate on vectors and matrices. Vectors are actually treated as one-dimensional matrices in Matlab. It is possible to use arrays of three and more dimensions, but we do not need them.
You can define an arbitrary row vector by typing
∗Course on Stochastic Processes for Finances 2016
1
x=[1,3,4]
Replace commas by semicolon and you get a column vector
y=[1;3;4]
You can also do it by transposing vector x
y1 = x’
Regularly spaced row vector t with values between [a,b] can be created with predetermined steps dt as
t = (a:dt:b)
or with predetermined number of steps N
t = linspace(a,b,N)
Note that the operation ended with semicolon ’;’ is not printed in Command Window, but those ended with comma or without punctuation will be printed. This is useful to see your results.
Element-wise multiplication between two matrices of the same size is done with operator ”.∗”.
C = A.∗B
Matrix multiplication of two matrices A, m × n, and B, n × k, is simply C = A∗B
and results in matrix C, m × k.
Exercise 1 (Cholesky decomposition). To practice a bit with matrices let us look at Cholesky decomposition, which will be used later to correlate random variables. Given a symmetric positive definite matrix C, Cholesky decomposi- tion factorizes it into a product of upper triangular matrix U and its transpose lower triangular matrix UT
Define correlation matrix C
C = UT U (1)
1 0.4 0.3 C = 0.4 1 0.5,
0.3 0.5 1
compute matrix U using Cholesky decomposition U = chol(C);
and show that Eq. 1 is fulfilled.
2
1.4 Graphs
Matlab has very powerful graphical capabilities. We will use them to visualize the results. To plot vector x as a function of its index, for example, one can type
plot (x)
A plot of vector x as a function of another vector t can be done as
plot(t ,x)
A plot of all rows in m×n matrix Y as a function of 1×n vector t is simply
plot (t ,Y)
It is often convenient to plot different graphs in different windows. Put figure(i)
in front of the plot command to direct the plot in figure i
figure (1) plot(t ,x) figure (2) plot (t ,Y)
The first plot will be window 1, while the second will be in window 2. Add legend if there are more than one line on the plot.
plot (x , sin (x) ,x , cos (x)) legend ( ’ Sin ’ , ’ Cos ’ )
In some exercises we need to plot multiple figures. It may be convenient to combine them in tabs of one figure. The following construction can be used.
It is a bit cumbersome but let us to avoid lots of standalone figures. For example in Exercise 9, you can write
hfig9 = figure (9); %Declare figure 9 and get handler ’hfig9 ’ set ( hfig9 , ’Name’ , ’ Geometrical Brownian Motion ’ ) ; %give name set ( hfig9 , ’ WindowStyle ’ , ’ normal ’ ) ; %make standalone window htabgroup9 = uitabgroup(hfig9 ); %create a group of tabs
htab9 1 = uitab ( htabgroup9 , ’ Title ’ , ’Abs returns ’ ) ; hax9 1 = axes(’Parent’, htab9 1);
plot(t,x) %plot x(t) in tab ’Abs returns’
htab9 2 = uitab(htabgroup9 , ’Title ’ , ’Log returns ’); hax9 2 = axes(’Parent’, htab9 2);
plot(t,y) %plot y(x) in tab ’Log returns
% etc.
These are just suggestions. Use your favorite method to visualize your results.
3
2 Random number generation 2.1 Random numbers in one dimension
We will use normally distributed random numbers a lot. The function to gen- erate a matrix n × m of independent random numbers is
Z = randn(n,m);
This function generate pseudo-random numbers of standard normal distribution, i.e. with mean zero and variance 1.
Exercise 2 (Random numbers). Generate vector, say Z, of normal random numbers of length N = 100, 000. Plot this series by using
figure (1) plot (Z)
Remember declaration ‘figure(1)’ above is not required, but is convenient in case you want to plot the next graph in your code in a different window. Put figure(2) in front of the second figure and so on. You better plot every 100-th value to make the plot ‘nicer’. You can do it by putting Z(1 : 100 : end), which means that only indexes from 1 till the last one with step 100 are used. Compute mean, standard deviation, excess kurtosis, and skewness of the series, see that these are indeed standard normal parameters:
mZ = mean(Z)
sZ = std(Z)
ekZ = kurtosis(Z) − 3 skZ = skewness(Z)
Observe that mean(Z) is close but not quite equal to zero due to the finite sample. The deviation will be even higher if we take a half of vector Z
ZH = Z(1:N/2) mean(ZH)
Combine ZH with itself with negative sign and compute the mean.
ZA = [ZH;−ZH];
This so called ‘antithetic sampling’ is a legitimate technique to reduce vari- ance in random number simulations with symmetric pdf. Why does this trick ‘magically’ make the mean exactly zero? Would this ‘magic’ work for any anti- symmetric component of the averaging function, like x3 for example?
Plot a histogram of the empirical density function. Use the following code
figure (2)
[y,x] = hist(Z,100); plot (x , y)
4
Function hist sorts elements in Z into 100 equally spaced bins. y is the number of counts in bins, while x are the center values in each bin. Function plot draws the graph of y as a function of x.
Modify your code to compare the obtained distribution with the standard nor- mal distribution, computed as
yn = normpdf(x); p l o t ( x , y , x , yn )
Why do the plots not coincide? How do you need to rescale y to match normal pdf yn?
In practice you need to generate vectors with non-zero mean and non-unit vari- ance. Use the very same vector Z, computed in the previous step, to obtain another vector Z1 with mean 3 and standard deviation 2. Compare to the corresponding normal pdf (do not forget to rescale y).
yn = normpdf(x,3 ,2); p l o t ( x , y , x , yn )
2.2 Random numbers in many dimensions
Usually, we need to generate random numbers for two or more correlated random variables. We use the fact that the sum of two normal random variables is again a normal random variable.
Suppose we have realizations of three independent random variables X, Y , and Z, all with zero mean. We can generate three correlated random variables, X ̃, Y ̃ , and Z ̃, by making linear combinations of vectors X, Y , and Z. Or in matrix form, we multiply matrix [XY Z], which consists of three columns X, Y , and Z,
x1
x2 [XYZ]= .
y1 z1 y2 z2
. . . …
xN yN zN
with some loading 3 × 3 matrix U :
X ̃ Y ̃ Z ̃ = [ X Y Z ] · U .
It turns out that in order to have a correlation matrix C between X ̃, Y ̃, and Z ̃, the corresponding matrix U is obtained by the Cholesky decomposition.
It can be shown analytically that the random variables produced in this way have indeed correlation matrix C. We will do it empirically.
For that you can use the standard formulas. Suppose that X = {x1, x2, …, xN } and Y = {y1 , y2 , …, yN } are series of observations for two normal random vari- ables with mean zero. We can estimate the covariance and correlation between
5
them in the following way
cov[X, Y ] = E[X · Y ] ≈ N − 1
or by using corresponding Matlab functions.
Exercise 3 (Independent normal random variables). Use randn(N, M ) to gen- erate a two column matrix Z with say N = 10, 000. Check that the columns are linearly independent. For that, compute covariances and correlations either by
cov Z = cov(Z(: ,1) ,Z(: ,2)); rho Z = corr(Z);
or
cov z = Z(: ,1) ’∗Z(: ,2)/(N−1);
rho z = Z(:,1)’∗Z(:,2)/sqrt((Z(:,1)’∗Z(:,1)) ∗ (Z(:,2)’∗Z(:,2)));
Think what exactly is done in the second method and why it gives slightly different values than the first method?
Plot Z(:, 1) vs. Z(:, 2) and observe a round ‘cloud’ of points.
Exercise 4 (Correlated normal random variables). Define in Matlab correlation matrix
1 −0.7 C=−0.7 1 ,
Compute matrix Z1 with two columns of correlated random variables using matrix Z and the Cholesky decomposition.
Z1 = Z ∗ U;
Compute covariance and correlation matrices, check that the correlation matrix is C. Plot Z(:, 1) vs. Z(:, 2) and observe that the cloud becomes elongated. Feel free to play with positive and negative correlations.
Exercise 5 (If Cholesky fails). Define correlation matrix
1 0.9 −0.9 C = 0.9 1 0.95 ,
−0.9 0.95 1
try to apply the Cholesky decomposition. The algorithm to compute the de- composition is very robust. If it fails it means that the matrix is not positive definite. Think in terms of correlations why the decomposition fails? Are all correlation matrices possible?
6
1 N i=1
xiyi
Ni=1 xiyi
E[X · Y ]
ρXY = E[X2]E[Y 2] ≈ Ni=1 x2i · Ni=1 yi2
3 Monte Carlo simulation of stochastic processes 3.1 Brownian motion generation
We start with generation of paths for Brownian motion W(t). Remember that
W(0) = 0
E[W(t) − W(s)] = 0
var [W (t) − W (s)] = t − s
E [(W (t) − W (s)) (W (v) − W (u))] = 0, ∀s ≤ t ≤ u ≤ v
We simulate Brownian motion on the interval [0,T]. For that we introduce a time grid {t0 = 0,t1 = ∆T,t2 = 2∆T,…,tM = M∆T = T}. At any time ti the value of the process for a single realization can be presented as
ii
W (ti) = (W (tj ) − W (tj−1)) = ∆W (tj ) j=1 j=1
Here ∆W(tj) is a normal random variable √
Zj∼N(0, ∆T)
√
∆T.
j=1 j=1
Exercise 6 (Single path generation). Generate and plot a single path of Brow- nian motion on [0, 1] interval. Try it a few times to see how it behaves from a realization to realization.
Exercise 7 (Multiple paths generation). Suppose we want to compute an ex- pectation of some function of Brownian motion at time T. We need to generate multiple paths. For that we need to loop over time steps and over different paths. Or we can use the Matlab matrix functionality to simplify the code and to speed up the calculations. There are different ways but the most efficient and transparent way is to generate a vector of increments for N paths on every time step and loop over the time steps like shown below.
W= zeros(npath,nstep+1); % Initial initialization of % the simulated variable (to zero in this case ).
for i = 1: nstep
Z = randn(npath ,1);
% Do what you need to do here end
WT=W(:,end); % Values at time T 7
∆T. Thus if {Zj} is a series of stan- ii
with mean zero and standard deviation
dard normal random numbers, we can generate a single realization for W(ti)
as
√ W(ti) = ∆W(tj) = Zj
Plot the paths to visualize them for yourself. You can do it in one go. If tgrid is a 1 × (nstep + 1) time grid vector, then plot
plot(tgrid , W)
Observe a square root divergence of the paths.
Generate N = 10, 000 paths and show that W (t = 1) has indeed zero mean, vari- ance one, and zero excess kurtosis and skewness. Plot the empirical distribution of W(T) at time T and compare it to the normal pdf.
Exercise 8 (dW2 = dt and dW1dW2 = 0). Show numerically that the following relations are satisfied
t 0
t
I1(t)= I2(t) =
dW2 =t dW1dW2 = 0
0
For this, approximate integrals I1 and I2 by sums, respectively,
i
S1(ti) = ∆W(tj)2 j=1
i
S2(ti) = ∆W1(tj )∆W2(tj ) j=1
Here ∆W1,2(tj) are independent normally distributed increments. Simulate 10 paths on the time interval [0,1] with different number of time steps N: 10, 100, . . . , 105. Plot the paths for each number of steps. See if the paths for Si converge to S1(ti) = ti and S2(ti) = 0, respectively.
3.2 Itˆo lemma
GBM is described by the following stochastic differential equation (SDE): dS(t) = rS(t)dt + σS(t)dWt
S(0) = S0 or by the discrete-time approximation
S(ti+1) = S(ti) + rS(ti)(ti+1 − ti) + σS(ti)(Wi+1 − Wi) Using Itˆo lemma the equation can be written as
σ2 dln(St)=(r− 2)dt+σdWt
or in the discrete-time form
s(ti+1) = lnS(ti+1) = lnS(ti) + (r − 2 )(ti+1 − ti) + σ(Wi+1 − Wi)
σ2 Observe the additional term σ2/2 in the drift.
8
Exercise 9 (Prove numerically Itˆo lemma). Simulate S(ti) and s(Ti) = ln S(ti).
1. Show that the distributions of S(tN ) and Sˆ = exp(s(tN )) are the same.
2. Compute log-returns lnX(tN)/X(t0) for both processes, S(ti) and Sˆ(ti) and show that they both have the same normal distribution.
You can use the following parameters for the model:
S0 = 1 r=0
σ = 0.2
You can also plot paths for S(ti) and exp(s(ti)) and observe that they are close but there is a difference. Simulation of ln S is exact in this case, while simulation of S(ti) converges to exp(s(ti)) only for small time steps.
3. Plot log-returns as a function of time for a single path. Take 1000 time steps to obtain a dense path. The volatility is pretty constant visually as one can expect.
3.3 Other stochastic processes
GBM model predicts log-normal returns as we could see in the previous exercise. This is typically not the case in practice. Look at the log-returns of DAX index and observe the volatility clustering, i.e. splitting in regions where volatility is small and regions where it is high. This indicates that volatility does not stay
constant in time. To account for that more sophisticated models are used. One class of them is called stochastic volatility models. There are different way to introduce stochastic volatility. Two very popular models are Heston and SABR. The first is popular in equity derivatives, the second in interest rate derivatives.
9
3.3.1 Heston model
Heston model has a mean reverting variance. The process for stock value S and variance ν are given by the following SDE.
dS(t) = rS(t)dt + ν(t)S(t)dW1(t)
dν(t) = κ(θ − ν(t))dt + ξν(t)dW2(t) ⟨dW1, dW2⟩ = ρdt
S(0) = S0 ν(0) = ν0
Here r is a risk-free rate, κ is a speed of mean reversion, θ is a long term variance, and ξ is a volatility of volatility. Brownian motions W1 and W2 are correlated with correlation ρ.
3.3.2 SABR model
The stochastic α, β, ρ (SABR) model is given by the following SDE
dF(t) = σ(t)F(t)βdW1(t) dσ(t) = ασ(t)dW2(t)
⟨dW1, dW2⟩ = ρdt F(0) = F0 σ(0) = σ0
Here F is a forward rate, usually defined in case of stock as F = F(t,T) = Ser(T−t). Forward rate is a rate at which two parties agree to buy/sell stock at a future time T. α is a volatility of volatility. β is some power of Forward rate that controls the skewness of return. Brownian motions are again correlated.
Exercise 10 (Stochastic volatility). Choose one of the models, Heston or SABR. Simulate the paths to maturity T = 1 year.
1. Plot a single path and look if you can observe a volatility clustering.
2. Plot histograms of log-returns and compare them to the normal distri- bution with the same mean and variance. Observe the skewness and fat tails.
10
As initial parameters for the model choose the following:
Heston SABR S0 = 1 F0 = 1 r=0 r=0
ν0 =0.22 θ = ν0
ρ = −0.7 ξ = 0.5 κ=1
4 Option pricing
σ0 =0.2 β = 1 ρ = −0.7 η = 0.5
In the current computer practical, we will use Monte Carlo simulation to price options on stock. We start with the celebrated Black-Scholes model and continue with other models. In all models, we will follow almost the same route to emphasize similarities and differences between models.
A small remark on the simulation approach. In the previous computer practical we simulated the stock price directly. We will continue to do this also in the current practical for clarity of exposition. However in many real life situations, it is better to use some transformations to make the computation quicker and more robust. See Appendix for performing computation in log-space, for example. The price of European type options, like call and put, can be computed using expectation under the risk-neutral measure as
V(0)=e−rTEQ0 [Φ(ST)].
Here St is a random process of stock price, ST is a stock price at maturity T , and Φ(ST ) is a payoff function. r is a risk-free rate. Let us express this expectation in terms of integral.
−rT Q −rT∞ V(0)=e E0 [Φ(ST)]=e
0
where φ(ST ) is a probability density function (PDF). In most cases, we do not know the analytical expression for PDF. We can, however, approximate it by an empirical PDF obtained by Monte Carlo simulations. MC simulation gives discrete points at maturity T, Si, and the corresponding empirical PDF can be expressed as a sum of Dirac delta functions δ(x − Si)
1 N
φe(ST ) = N
δ(ST − Si).
Φ(ST)φ(ST)dST,
i=1
11
The empirical PDF converges to the actual PDF as N goes to infinity in the sense that the following integral converges to the desired expectation.
∞ ∞
Φ(ST )φe(ST )dST → Φ(ST )φ(ST )dST .
N→∞ 0 Thus, we can approximate the expectation as
0
V(0)=e =e
=e N Φ(Si), i=1
−rT Q −rT∞ E0 [Φ(ST)]≈e
− r T ∞ 0
1 N
Φ(ST)N − r T 1 N
i=1
Φ(ST)φe(ST)dST δ(ST −Si)dST
0
i.e. by computing a discounted mean of the payoff function with the generated values of stock price Si.
4.1 Black-Scholes model
4.1.1 Call and put prices, comparison to Black-Scholes formula
Exercise 11 (Black-Scholes). Compute European call and put prices using Monte Carlo simulation code from the previous practical for Geometrical Brow- nian motion under the risk-neutral measure
dSt = (r − c)Stdt + σStdW S(0) = S0,
where r is a risk-free rate, c is a continues dividend rate, σ is a volatility. The
payoff functions are
(S − K)+ call, Φ(S) = (K − S)+ put.
The expectation is approximates as shown above
i=1
The first thing one usually checks, if the code produces the right results for simple cases. Compute the forward price of the stock F(T) at maturity T = 1 and compare it to the analytical expression
F (T ) = e(r−c)T S(0).
12
V (0) = e N
Φ(Si).
− r T 1 N
4.1.2 Put-Call parity
Graph 1: Check the put-call parity. Plot the graph for the difference between call and put prices C(K) − P(K) for different strikes K and compare it to e−cT S(0) − e−rT K.
4.1.3 Comparison to Black-Scholes formula
In case of GBM, we can use analytical solution to compare to our Monte Carlo simulations. The celebrated Black-Scholes formula
V (0) = ψ e−cT SN (ψd1) − e−rT KN (ψd2)
ln S/K + (r − c)T + σ2T d1= √ 2
√σT d2=d1−σ T
+1 for call −1 for put.
is given in the function bsprice(S,K,r ,c,v,T,pc)
The function takes as input the current stock price S, a vector of strikes K, risk-free rate r, dividend rate c, volatility v, time to maturity T, and put or call indicator (-1 for put and +1 for call).
Graph 2: Plot on the graph the Monte Carlo prices of call and put for different strikes together with the prices of Black-Scholes formula.
4.1.4 Implied volatility
It is market standard to quote the price of options in terms of the volatility. For this, traders use the Black-Scholes formula. All other parameters are usually readily available on the market. Moreover, the stock price is most volatile, and by quoting volatility to the counterparty traders avoid the dependency on the stock price. When the deal is finalized, the actual stock price at that moment is put into the formula.
Graph 3: Implied volatility. We will use the Black-Scholes formula to compute the volatility implied from your MC call and put prices and plot it as a function of strike. Function
bsimpvol (S ,K, r , c ,T, pc , price )
does this. The function takes the same arguments as ’bsprice’, except for the volatility, plus the prices of the options and returns the corresponding implied volatility. Compare on the graph the implied volatility and the Black-Scholes input volatility that we used in our simulations for different strikes.
13
ψ=
4.2 SABR model
Use the Monte Carlo framework for SABR developed in the first practical to compute option prices.
Exercise 12 (SABR). 4.2.1 Put-Call parity
Check that the simulated forward price is in line with the analytical formula. Graph 1: Do as in the Black-Scholes case a check that the put-call parity is satisfied.
4.2.2 Call and put prices, comparison to Black-Scholes formula
Graph 2: Compare the simulated prices to the Black-Scholes formula. Use the volatility of log-returns in Black-Scholes formula. You will see the price difference and slightly different shape. It is not convenient to compare prices as they depend strongly on strike. We will use implied volatility for more proper comparison.
4.2.3 Implied volatility
Graph 3: Compute the implied volatility for SABR model. Plot it together with the Black-Scholes volatility as a function of strike. Observe the smile and the skew of the volatility σimp(K). Play with different parameters of the SABR model to observe which parameters are responsible for the level, slope (skew), and curvature (smile) of the volatility smile.
4.2.4 Volatility approximation
The reason why the SABR model is so popular is the existence of an analytical approximation for implied volatility in SABR model. This makes the calibration of the model to the market quotes very easy. The derivation of this approxima- tion is very technical [1], and is beyond this course.
Graph 4: Use function ’sabrVolatility(S,K,r,c,alpha,beta,rho,eta,T)’ to compute the implied volatility using this approximation and compare to the implied volatilities obtained from your simulations.
4.3 Heston model
Use the Monte Carlo framework for SABR developed in the first practical to compute option prices.
In case of the Heston it is not that easy to guaranty the positivity of the volatil- ity. There are different computational schemes to do it, which are beyond this
14
course. We will use a simple truncation to achieve this goal, i.e. we will use v+ = max(v, 0) in our simulations.
dSt = (r − c)Stdt + νt+StdW1 dνt = κ(θ − νt+)dt + ξνt+dW2
⟨dW1, dW2⟩ = ρdt
where we truncated variance ν to positive values ν+. Exercise 13 (Heston).
4.3.1 Put-Call parity
S(0) = S0 ν(0) = ν0
Check that the simulated forward price is in line with the analytical formula. Graph 1: Do as in the Black-Scholes case a check that the put-call parity is satisfied.
4.3.2 Call and put prices, comparison to Black-Scholes formula
Graph 2: Compare the simulated prices to the Black-Scholes formula. Use the volatility of log-returns in Black-Scholes formula.
4.3.3 Implied volatility
Graph 3: Compute the implied volatility for the Heston model. Plot it together with the Black-Scholes volatility as a function of strike. Observe the smile and the skew of the volatility σimp(K). Play with different parameters of the Heston model to observe which parameters are responsible for the level, slope (skew), and the curvature (smile) of the volatility smile.
4.4 Local volatility model
SABR and Heston are parametric models and in many cases do not allow to fit all market quotes for all strikes and maturities exactly. There is another type of model, called Local Volatility model, which allows to fit all the available quotes on the market . In this model the volatility is a deterministic function of stock price and time, σloc(St, t):
dSt = (r − c)Stdt + σ(St, t)StdW S(0) = S0.
Generally, it is a daunting task to compute the local volatility function for real quotes. We however will follow an easy way and assume that the market quotes
are completely in line with the SABR model. Then we can use the SABR ap- proximation for the Black-Scholes implied volatility to compute the local volatil-
ity function. The corresponding function ’sabr2LocalVol(S0,K,r,c,alpha,beta,rho,nu,t)’ returns the local volatility σ(St = K, t) for a given vector of strikes K and ma-
turity t.
15
Again we perform the Monte Carlo simulation is as in case of GBM. Only instead of the constant volatility, we put in the calculations the local volatility, corresponding to the current time step and the current stock price at this time and in this path.
Exercise 14 (Local volatility).
4.4.1 Call and put prices, comparison to Black-Scholes formula
Graph 1: Calculate call and put prices and compare them to the Black-Scholes. Plot also the call and put prices computed with the volatility given by the SABR approximation.
4.4.2 Implied volatility
Graph 2: Compute the implied volatility and compare it to Black-Scholes and to the SABR approximation.
4.5 ⋄ Pricing Barrier option
To show how different models affect pricing of exotic derivatives we look at barrier option pricing in SABR and local volatility models. These models give the very same European call and put prices for all strikes and maturitites. So, we may expect that the barrier prices are the same as well, but they are not. This shows that the choice of the model is very important for pricing exotic options.
A barrier option has the same payoff as vanilla option, however an Out option can be knocked out if the stock price crosses a predefined barrier during the life time of the option. If the stock price crosses the barrier at any time during the life time, the option ceases to exist, and the holder does not get anything even if the vanilla option is in the money. The payoff of the up-and-out option for example is
Ψ(ST)=(ST −K)+I(St