MTHM003 ANALYSIS AND COMPUTATION FOR FINANCE
Coursework 2: Due Monday 28 November 2022
Overview and Submission information
General guidelines. Please write down your answers clearly quoting numerical values to 4 signif-
icant figures, unless otherwise stated. Full marks achieved in this assignment will contribute 15%
of the final module mark. Usual university penalties apply for late submission unless mitigation
is approved. Feedback will be written on your work, and general comments will also appear on
ELE.
Submission format. Please submit a single PDF file summarising your answers to the questions
(e.g., computed quantities, graphs, plots, etc.) and including all supporting MATLAB code via
eBART submission. Any queries, please ask at the front desk for help. Analytical calculations
and discussions (when asked for) can be done by hand-writing.
Allocation of marks. When solving any particular problem, marks will be awarded for clearly
commented (MATLAB) codes, and carefully reasoned mathematical solutions. Any diagrams,
graphs, etc., should be clear to understand, with good comment or clear labelling. For questions
requiring numerical simulation/output, marks will primarily be awarded for demonstrating use of
MATLAB: little credit will be given for other solution schemes.
Academic Misconduct. All work submitted including functions and scripts must be your
own and not copied from other students or other sources (e.g., web pages). Any evidence
of possible plagiarism will be investigated following university procedures. You are welcome to
discuss ideas and possible methods, of course, with each other and with lecturers, etc., but the
work you hand in must be your own.
Main contact for this assignment: Frank Kwasniok (F.Kwasniok@exeter.ac.uk)
Questions
1. Chaotic dynamics in a simple system of ODEs.
(a) (10 marks) Write a MATLAB function RK4.m to solve a system of ODEs using the 4th
order Runge–Kutta method. Recall that to obtain the solution at time tj+1 from the state
yj at time tj, this iterative scheme requires you to calculate four variables
k1 = f(yj, tj),
k4 = f(yj + hk3, tj + h),
and time step using
yj+1 = yj +
h
6
(k1 + 2k2 + 2k3 + k4).
Your implementation should use the same inputs and outputs as the explicit Euler and
Euler–Heun methods given in the lecture notes.
(b) (5 marks) The Ro¨ssler equations are a system of ODEs known to produce chaotic dynamics.
The equations are related to the Lorenz equations describing atmospheric convection. For
a column state vector X = (x, y, z)T, time t and parameters a, b and c the equations are
defined as follows:
x˙ = ?y ? z,
y˙ = x+ ay,
z˙ = b+ z(x? c).
For c2 > 4ab the equations have two fixed points X± = (x±, y±, z±)
T given by
x± = ?ay±,
y± = ?z±,
z± =
c±√c2 ? 4ab
2a
.
(i) Fix the parameters to a = ?0.1, b = 0.2 and c = 5.7. Starting from the initial
condition X? = (0, 0, 0.25)T numerically solve the Ro¨ssler equations using your RK4
implementation for a time interval t ∈ [0, 200] with a step size of h = 0.005 (use the
same step size for all questions). Plot each of the state variables against time.
(ii) You should find that the trajectory converges to one of the fixed points, give the value
of this fixed point to 4 decimal places as calculated in MATLAB through simulation;
check your computed value against the above formulae.
(c) (10 marks)
(i) For the parameter values a = 0.1, b = 0.2 and c = 5.7, the same initial condition as in
(b), and a time interval t ∈ [0, 800], numerically solve the Ro¨ssler equations. Discard
any transient behaviour (say t < 300) and plot the trajectory in (x, y, z) coordinates.
You should find a closed curve (a periodic orbit).
(ii) Taking the Poincare′ section P = {x, z : y = 0, x > 0} (the positive-x part of the plane
y = 0), find all points at which the trajectory passes through P . Plot the trajectory in
(x, z) coordinates, marking any points along the trajectory that cross through P . At
a = 0.1 these points should all coincide exactly. Give the value of the periodic orbit’s
minimal period T to 4 decimal places.
(iii) Plot each of the state variables against time, marking the points where the trajectory
crosses through P .
(d) (5 marks) Repeat part (c) for a = 0.13, a = 0.15, a = 0.1525 and a = 0.2. Comment on
how the dynamics change as a increases.
Hints:
1. The Ro¨ssler equations can be defined as an anonymous function in MATLAB to work with
your function RK4.m as follows:
ross=@(t,x)[...] % function defining right-hand side
2. You can check the equations are defined correctly by calling ross(0,Xpl) or ross(0,Xmi)
with Xpl and Xmi being the fixed points given above, which should return the column vector
(0, 0, 0)T.
3. In part (c) the points should coincide exactly. You can improve the accuracy by interpolating
from points along the trajectory either side of P to find more accurate values of x and z
when crossing through at y = 0.
[30]
2. Fokker–Planck equation.
A simple model for asset prices St for t ≥ 0 is given by the following stochastic differential equation
(SDE):
dSt = μStdt+ σStdWt, S0 = s0,
where Wt is a standard Wiener process, μ is the drift coefficient or expected return and σ is
the volatility. You are provided with a MATLAB script sde example cw2.m that generates 500
realisations of this SDE and plots the distribution of the trajectories at final time T = 2 as a
histogram for the parameter values μ = 0.2, σ = 0.25 and initial conditions s0 drawn from a
normal distribution N(ν, λ2) centered at ν = 0.5 with standard deviation λ = 0.05.
The probability density p(t, s) of trajectories of the above SDE over the spatial variable asset price
s ∈ (0,∞) satisfies the Fokker–Planck equation (FPE)
tp(t, s) = ??s[μsp(t, s)] + 1
2
ss[σ
2s2p(t, s)],
a partial differential equation. The FPE can be solved numerically using a forward in time, central
in space (FTCS) finite difference method as described below:
Discretise the time interval t ∈ [0, 2] as tj = t0 + j?t for j = 0, 1, . . . , nt, where the final
time point is T = t0 + nt?t = 2.
For s ∈ [0, 3] let sk = (k ? 1)?s with k = 1, 2, . . . , ns where ns = 3/?s + 1. While s can
take any value in (0,∞), here we truncate the solution assuming p(t, s) = 0 for s ≥ 3 for
the given time interval.
Notation: pjk ≈ p(tj, sk) for j = 0, 1, . . . , nt and k = 1, 2, . . . , ns
Initial condition (IC): p0k = p(0, sk) = (1/
√
2piλ) exp[?(sk ? ν)2/(2λ2)] for k = 1, 2, . . . , ns,
that is, the normal distribution N(ν, λ2)
Dirichlet boundary conditions (BCs): pj1 = p(tj, s1) = 0 and pjns = p(tj, sns) = 0 for
j = 0, 1, . . . , nt
Recursively define (loop over j and k):
p(j+1)k = pjk μt
2s
[sk+1pj(k+1) sk1pj(k1)] + σ
2t
22s
[s2k+1pj(k+1) 2s2kpjk + s2k1pj(k1)]
(a) (20 marks) Write MATLAB code for the above FTCS method and simulate the FPE with
t = 0.0001, ?s = 0.01, parameter values, initial condition and boundary conditions as
given above.
(b) (10 marks) Demonstrate that your implementation works and is accurate by producing the
following plots:
– Plot your solution p(t, s) as a surface or as a 2D colour map (use, for example, surf or
imagesc).
– Plot the distributions p(t, s) at the initial and final time points, and compare these
with the distribution of St at the final time T = 2 as produced by the given script
sde example cw2.m.
Hints:
Your solution to the FPE will be a matrix p of size (nt+1)×ns. Before using the recursive equation
given above in a for-loop you can pre-allocate some values in p using IC for the first row, and
using BCs for the first and final columns. [30]
3. The diffusion equation and the Black–Scholes equation.
(a) (15 marks) Use the MATLAB PDE solver pdepe to solve the diffusion equation
on the interval x ∈ [0, 1] with the initial condition
u(0, x) = u0(x) = ? sin 5pix
2
subject to the boundary conditions
u(t, 0) = 0
and
(t, 1) = 0.
Use 41 equally spaced mesh points in the x-direction and integrate between t = 0 and t = 0.1
using 21 equally spaced time steps, and follow the method given in the example in section 12
of the lecture notes (edit the functions pdefun, pdeic and pdebc for use with the MATLAB
function pdepe). For output, plot u(t, x) at the times t = 0, t = 0.025, t = 0.05, t = 0.075
and t = 0.1 on one graph.
(b) (25 marks)
(i) Solve the transformed Black–Scholes equation in the European put case
with initial condition
u(0, x) = u0(x) = max
(ii) Solve the transformed equations for the case E = 10, r = 0.1 and σ = 0.45 with
T = 1/3, where k = 2r/σ2. Explain why the scaled time τ runs from 0 to 0.03375.
Your output should be P at integer values of S from S = 0 to S = 16 (noting that
P0 = Ee
kτ can be determined analytically; explain why this is the case). Some of
these values are given in the table below for comparison.
S P
0.00 9.6722
2.00 7.6722
4.00 5.6723
6.00 3.6977
(iii) Experiment by increasing the number of mesh points in the x-direction, and by in-
creasing the value of xb, which gives the points where the boundary conditions are
applied.
It is important to annotate your code with comments, explaining what each of the function
m-files does and how these function m-files refer to the general equation solved by pdepe.
Annotate also the driving m-file script, explaining in particular the role of xb and the use of
pdeval in extracting P in terms of S for the output.