I. Background
The purpose of this practical is to investigate by simulation approximations used in stochastic
integrals and SDEs. The R may be found at All Programs / Statistics ... The following functions
in R may help to address some of the questions. These can be read by clicking on File ! Source
on the R menu bar, at Prac.R, which opens another R window where the le may be then copied
and pasted into the R working window; the address of the le is (computer/ network locations /
n:/lucs-h13/userlib/userlib/stalib3/math3733/prac.R Alternatively, you may use the copy of this
le found on the VLE at Learning resources/Practical. Remember that you must read { by
copying it in your R Console window { the Prac.R le BEFORE you try examples or
start working on your main task. Also, you yourself must read all lines of the "Important
remarks" section below before trying with R.
PLEASE, DO NOT WORK IN THE DIRECTORY n:/lucs-h13/userlib/userlib/stalib3/math3733/
Some important remarks (to read very carefully in the very beginning!)
1. bmsim(nsim,ntime) This function creates an nsim (ntime + 1) matrix. Each row is an inde-
pendent simulation of a standard Brownian motion, observed at the time points t = j=ntime; j =
0;:::;ntime. To avoid memory problems, it is recommended that you limit the dimensions so that1
nsim 1000; ntime 1000.
2. stoch.int(fun1=fun.null, fun2=fun.null,w=bm, lh=T, sigma=1, nu=1, alpha=1) This
function computes an approximation to the combination of stochastic & Riemann integrals (writing
f1 and f2 for fun1 and fun2)
Z 1
0
f1(Wt)dWt +
Z 1
0
f2(Wt)dt
(apparently, with = = 1 in this example). If lh=T (i.e. left-hand = \true", the default option),
then the functions f1 and f2 are evaluated at the lefthand endpoint of each interval [j=n; (j+1)=n],
and the approximation used is (writing n for ntime)
n 1X
j=0
f1(Wj=n) (W(j+1)=n Wj=n) + 1n f2(Wj=n)
;
whereas if lh=F (i.e. left-hand="false") the approximation becomes
n 1X
j=0
f1(W(j+1)=n) (W(j+1)=n Wj=n) + 1n f2(W(j+1)=n)
1Depending on the computer, 5000 might be tried.
MATH 3733 Practical Page 2
The approximation is evaluated for each of the nsim simulations, and the answer is output as a
vector of length nsim.
3. sde(fun1=fun.null, fun2=fun.null,w=bm, init=0, sigma=1, nu=1, alpha=1) This func-
tion constructs a discrete approximation to the solution Xt of the stochastic di erential equation
dXt = f1(Xt)dWt + f2(Xt)dt
using the method decribed in the lectures. The result is output as a matrix of size nsim ntime+1,
the same dimensions as the input Brownian motion. Each row contains an independent simulation,
giving the values of the process at the times t = j=n; j = 0;:::;n.
4. Several functions have been provided to use in stoch.int and sde. Each is a function f(x) of
a real-valued argument x. Take care when you choose a proper function for your work: it is your
task to decide, which function is needed.
(a) fun.null, the null function, f(x) = 0.
(b) fun.one, the constant function, f(x) = 1.
(c) fun.id, the identity function, f(x) = x.
(d) fun.sq, the square function, f(x) = x2.
(e) fun.exp, the exponential function, f(x) = exp( x). The choice of is input as an argument
to stoch.int or sde.
4. dw2(w=bm) This computes the sum of the squared increments Pn 1j=0 (W(j+1)=n Wj=n)2 for each
of the simulations and outputs a vector of length nsim.
5. There are standard functions in R to study the properties of a vector of data, given in y, say.
(a) hist(y), plots a histogram
(b) qqnorm(y), constructs a qq-plot of the data quantiles vs. the population quantiles for the
normal distribution. If the data are normally distributed N( ; 2), the plot should be ap-
proximately a straight line passing through (0; ) and of slope .
(c) qqchisq(y, df), a specially written function which constructs a qq-plot for a 2df distribution.
If the data follow a multiple 2 of this distribution, the plot should be a straight line of slope
2 passing through the origin. You are not supposed to perform. any hypothesis testing.
(d) mean(y), var(y), mean and variance of y; pnorm(a), gives Laplace function value (a).
(f) log(a), gives the value of lna; abs(a) gives jaj.
II. Examples
The following commands illustrate how to simulate a Brownian motion, to evaluate
Z 1
0
7 exp(2Wt)dWt +W2t dt
MATH 3733 Practical Page 3
and how to solve dXt = X2t dWt, starting at X0 = 6. Note that the latter equation does not satisfy any
existence/uniqueness conditions, which may result (or may not) in \explosion" of solution over a nite time
interval. If this occur, the programme may report an error or a problem.
bm<-bmsim(200,100) # create 200 simulations of BM, each at 100
times on [0,1]
a<-dw2() # sum the squared increments of BM for each simulation
mean(a) # print out the mean of a hist(a) # plot a histogram of a
si1<-stoch.int(fun1=fun.exp, fun2=fun.sq, sigma=7, alpha=2)
si2<-stoch.int(fun1=fun.exp, fun2=fun.sq, lh=F, sigma=7, alpha=2)
# two evaluations of stochastic integral with different
# endpoint conventions
hist(si1+si2) # histogram of si1+si2
x<-sde(fun1=fun.sq, init=6) # simulate solution of sde
x1<-x[,101] # vector of simulated values of X at time t=1
qqnorm(x1) # assess normality of x1
qqnorm(rnorm(100)) # simulate 100 N(0,1) variates and see what the
# qq plot looks like
III. The main task
1. Simulate a Brownian motion (BM) with nsim=1000, ntime=1000 and with a few other values
(not too big), store the result in a matrix bm1, say. Note that simulations for W1 are stored in
bm1[,ntime+1].
(a) Simulate a BM again and store it, e.g., in a matrix bm2. Now de ne bm3 = (bm1 + bm2)=p2.
Verify whether or not bm3[,ntime+1] is normal. What are its parameters mean and variance?
(b) How would you test whether bm3 is a standard Wiener process? (You are not asked to perform
this testing, just report in wording how would you do it.)
(c) Investigate the approximations to R10 2WtdWt using lh=T and lh=F; store the results in two
vectors, si1 and si2, say; for several values of ntime, investigate the following issues:
(1) how does the sample distribution of si1 compare to the theoretical distribution of the stochas-
tic integral? (It is your task to decide what is the theoretical distribution)
(2) how would you expect the di erence si1 - si2 to behave as ntime tends to in nity? And
how does the di erence behave in practice? Report whatever you observe.
2. For several values of ntime simulate an approximate solution to the SDE
dYt = YtdWt +Ytdt; Y0 = 1:
Compare the sample distribution of Y1 with the theoretical distribution. It is your task to decide
what is the theoretical distribution; yet, the hint is, rst of all, to solve this SDE theoretically.
3. For the BMS stochastic market St = S0 exp ( Wt +t) with S0 = 1 compute the price C0 of
the call option with T = 1 and K = 1:1 at t = 0 assuming interest rate r = 1. You may use the
MATH 3733 Practical Page 4
formula exp( 1)~ECT and simulate the expectation2 by several samples of (WT) and (CT), that is,
using the Law of Large Numbers (LLN), replace exp( 1)~ECT by the sample mean
exp( 1) 1N
NX
j=1
CjT:
Do not forget that expectation should be with respect to the implied probability measure.
Compare the result with the one obtained from the theoretical formula for C0.
IV. Deadlines and comments
The supervised sessions will be on Monday 12.03.2018, 12:00-14:00 at Chemical and Process Engineer-
ing Cluster (GR.06), and on Monday 12.03.2018, 15:00-17:00 at Social Sciences Cluster (10.11), and
you can come to one of these sessions depending on your timetable. You are very strongly advised to
start working individually well before the session, which is NOT a lecture, nor a workshop on R,
but a class where a lecturer and an assistant may answer some of your basic questions related to the
practical. The questions like \why this code does not work", or \where is an error in my code" are
NOT considered as basic ones. The practice shows that usually the code does not work if you did not
read carefully the introduction to the task.
At http://infoscreen.leeds.ac.uk/clustercode/index.php you may nd the reminder about the code
for computer clusters access outside normal times.
The deadline for submission of your reports is 14:00 Friday 23.03.2018 as a printed stapled copy (partly
handwritten allowed, although not advised) in the grey mailbox level 8 named B. Thorpe. (Note
that I am not a grader and I will be away during the Easter break, so, please, DO NOT submit to my
mailbox.) Late submission { except special circumstances that must be formally reported to the
Taught O ce { will result in a penalty (the School policy is -5% per any missed day). It is advised
to email to me and the grader in case of each late submission.
Your reports are a part of your exam and will NOT be returned. Please, do not ask about it.
If you wish so, you may see your reports graded and with the grader’s comments later in my o ce.
You may work together, but each report must be written independently3. If plagiarism is sus-
pected, then it may have serious consequences: see http://www.lts.leeds.ac.uk/plagiarism/
Do not forget to SIGN your reports and put your name, your student ID number and the
code of the module. The pages must be stapled and, please, do not use plastic folders. If
all the pages are not stapled properly, there is no guarantee that some of them may be missing: also,
non-stapled pages will be in all cases removed from a plastic folder after which there is no guarantee
that they do not separate. If it happens, the grader would not be responsible for this! Also, better
you put your name and your student ID on each page of your report.
All available marks will be eventually put on the VLE in the Gradebook of the module.
R is available on the Web. If you wish so, you may do your practical at home as long as you can run
R on your computer (see the rst page). For help about R there is a "Help" in the R worked place;
also see https://www1.maths.leeds.ac.uk/school/students/computing.html
Your report is expected to have about 6 { 12 pages, but if necessary a few more pages should be OK if
necessary for graphs (which are very welcome) and the most important parts of the code. There might
be a penalty, in particular, if the grader is in doubt about what you did and he cannot check your
code. Yet, you are not assumed to copy the whole auxiliary le with the hints of the code provided.
2Do not forget that you have to use implied probabilities.
3In particular, your data should not be identical by all means!