数学计算代写: Math 3733 Practical

数值科学计算

I. Background

Math 3733 Practical – spring 2018

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 file may be then copied and pasted into the R working window; the address of the file is (computer/ network locations / n:/lucs-h13/userlib/userlib/stalib3/math3733/prac.R Alternatively, you may use the copy of this file found on the VLE at Learning resources/Practical. Remember that you must read – by copying it in your R Console window – the Prac.R file 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

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)

􏰟1 􏰟1
σ f1(Wt) dWt + ν f2(Wt) dt

00

(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−1􏰚 1􏰛 􏰞 σf1(Wj/n)(W(j+1)/n −Wj/n)+ nνf2(Wj/n) ,

j=0
whereas if lh=F (i.e. left-hand=”false”) the approximation becomes

n−1􏰚 1􏰛 􏰞 σf1(W(j+1)/n)(W(j+1)/n −Wj/n)+ nνf2(W(j+1)/n) .

j=0
1Depending on the computer, 5000 might be tried.

MATH 3733 Practical Page 1

(to read very carefully in the very beginning!)

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 differential 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.

  1. (a)  fun.null, the null function, f(x) = 0.
  2. (b)  fun.one, the constant function, f(x) = 1.
  3. (c)  fun.id, the identity function, f(x) = x.
  4. (d)  fun.sq, the square function, f(x) = x2.
  5. (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 􏰜n−1(W − W )2 for each j=0 (j+1)/n j/n

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.

  1. (a)  hist(y), plots a histogram
  2. (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 σ.
  3. (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.
  4. (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 |a|.

II. Examples

The following commands illustrate how to simulate a Brownian motion, to evaluate

􏰟1
7 􏰘exp(2 Wt) dWt + Wt2 dt􏰙

0

and how to solve dXt = Xt2 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 finite 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 define bm3 = (bm1 + bm2)/ 2. 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 􏰝 1 2Wt dWt using lh=T and lh=F; store the results in two 0

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?

(2) how would you expect the difference si1 – si2 to behave as ntime tends to infinity? And how does the difference behave in practice?

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, first 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 3

√

(It is your task to decide what is the theoretical distribution)

Report whatever you observe.

MATH 3733 Practical Page 4 formula exp(−1)E ̃CT and simulate the expectation2 by several samples of (WT ) and (CT ), that is,

using the Law of Large Numbers (LLN), replace exp(−1)E ̃CT by the sample mean 1N

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 find 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 Office – will result in a penalty (the School policy is -5% per any missed day). It is advised to email to me and the grader b.thorpe@leeds.ac.uk 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 office.
  • 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 first 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 file 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!

exp(−1) 􏰞 CTj . N j=1

Related Posts