# Efficient GeometricBrownianMotionProcess Monte Carlo simulation

Following the answers in this post, I'm trying to implement something similar. If the GBM stays inside the corridor [L, U] between predefined times it should return 1 otherwise 0. I wrote the following code:

MC[X_, t_, Ti_, per_, L_, U_, r_, s_, n_, prob_] :=
Module[{data = 0, data1, data2, value, j},
data =
RandomFunction[GeometricBrownianMotionProcess[r, s, X],
{0, Ti[[Length[Ti]]] + per, .01}, n];
value = 0;
For[j = 1, j <= n, j++,
data1 = Table[Select[#, #[] >= Ti[[i]] && #[] <= (Ti[[i]] + per) &] & /@
{data["Paths"][[j]]}, {i, 1, Length[Ti]}];
data2 = Flatten[data1, 2];
value +=
If[Max[Table[data2[[i]][], {i, 1, Length[data2]}]] <= U &&
Min[Table[data2[[i]][], {i, 1, Length[data2]}]] >= L,
1, 0]];
N[1/n*value]];

MC[100, 0, {1, 2}, 1, 80, 120, 0.01, 0.15, 100, prob]


If the GBM stays inside [80,120] between the times [1,2] and [2,3], value should be 1 otherwise 0.

Now this code works, but it is really slow. For instance for 10,000 simulations it takes about 10 minutes. Is there a way to improve this?

I also want to calculate the probability that the GBM breaches the corridor between two discrete simulations. For example, if $\mathbb{P}[X_t,X_{t+1},L,U]=prob[\cdots]$ and $\mathbb{P}>0.5$, I want to set value to 0. In short, I'm looking for a way to apply a certain function prob to every simulated value inside my corridor.

Is there a way to do this without a loop? Can someone can give me a hint on how to do this in an efficient way?

thank you for your answer. Yes I'm trying to price a corridor option and the function prob I mentioned should account for the penetration probability. But I want to get the value of the option if there are more corridor periods than one.

To visualize what I mean consider the slightly altered plot you provided:

Module[{S0 = 100, r = 0.01, σ = 0.15, L = 80, U = 120,
T = 1, Δt = 0.01, nTimes = 100, nPaths = 10,
seed = 123, drift, diff, times, paths, containedPaths, PathToLine},
SeedRandom[seed];
times = Range[0, T, Δt];
PathToLine[path_] := Line[Transpose[{times, path}]];

{drift, diff} = {Exp[(r - σ^2/2) Δt],
diff = σ Sqrt[Δt]};

paths = Table[GBMPathCompiled[S0, drift, diff, nTimes], {nPaths}];
containedPaths = Select[paths, (And[Max[#] < U, Min[#] > L]) &];
Graphics[{PathToLine /@ paths,  {Red, Line[{{0.2, L}, {0.4, L}}],
Line[{{0.8, L}, {1, L}}], Line[{{0.2, U}, {0.4, U}}],
Line[{{0.8, U}, {1, U}}]}}, Axes -> Automatic,
AspectRatio -> 0.25]]


So the corridor should be active just between certain times $[T_i,T_i+per]$, where $T_i$ are the starting times and $per$ the length of the period. For the example above, this would be Ti={0.2,0.8} and per=

• I am somewhat puzzled over why this question has attracted three downvotes. It is intelligible, clear, and contains complete working code. Although the latter may be far from perfect, performance tuning questions generally have not been received so poorly in the past. Perhaps the downvoters would care to offer some constructive criticism for the OP. Mar 3 '13 at 18:44
• I've moved your "answer" to be an edit of the question. I know that it takes a little while to get used to the Q&A approach of StackExchange sites, but please understand that posts here aren't over discussion threads. Answer posts should only be used for actual answers. A belated welcome to Mathematica.SE nonetheless. Mar 7 '13 at 4:09

Here is a faster way. I've only coded a portion of your problem to provide the flavor.

GBMPathCompiled =
Compile[{{S0, _Real}, {drift, _Real}, {diff, _Real}, {nSteps, _Integer}},
FoldList[(#1 drift Exp[ diff #2]) &, S0,
RandomVariate[NormalDistribution[0, 1], nSteps]]];


Let's first visualize the situation

Module[{S0 = 100, r = 0.01, \[Sigma] = 0.15, L = 80, U = 120,
T = 1, \[CapitalDelta]t = 0.01, nTimes = 100, nPaths = 10,
seed = 123, drift, diff, times, paths, containedPaths, PathToLine},
SeedRandom[seed];
times = Range[0, T, \[CapitalDelta]t];
PathToLine[path_] := Line[Transpose[{times, path}]];

{drift, diff} = {Exp[(r - \[Sigma]^2/2) \[CapitalDelta]t],
diff = \[Sigma] Sqrt[\[CapitalDelta]t]};

paths = Table[GBMPathCompiled[S0, drift, diff, nTimes], {nPaths}];
containedPaths = Select[paths, (And[Max[#] < U, Min[#] > L]) &];
Graphics[{PathToLine /@ paths, {Blue,
PathToLine /@ containedPaths}, {Red, Line[{{0, L}, {T, L}}],
Line[{{0, U}, {T, U}}]}}, Axes -> Automatic,
AspectRatio -> 0.25]]


The double barrier penetration probability is obtained with

Timing[Module[{S0 = 100, r = 0.01, \[Sigma] = 0.15, L = 80, U = 120,
T = 1, \[CapitalDelta]t = 0.01, nTimes = 100, nPaths = 10000,
seed = 123, drift, diff, times, path, results},
SeedRandom[seed];
{drift, diff} = {Exp[(r - \[Sigma]^2/2) \[CapitalDelta]t],
diff = \[Sigma] Sqrt[\[CapitalDelta]t]};
results = Table[path = GBMPathCompiled[S0, drift, diff, nTimes];
If[ (And[Max[#] < U, Min[#] > L]) &[path], 0, 1.], {nPaths}];
{Mean[results], StandardDeviation[results]/Sqrt[nPaths]}]]

(* {0.171601, {0.322, 0.00467266}} *)


With a bit of work, this code could almost certainly be speeded up. Also, this can likely be done much more elegantly using the new MMA9 stochastic process capabilities.

Note: Using Max and Min on the path underestimates the penetration probability because you miss the penetrations that occur between sample times. You can fix this with a Brownian bridge.

Note: Of course, you are probably working up to a more complicated process. This double barrier problem with pure GBM has an analytic answer in the form of a rapidly converging infinite sum.

Then, you must build the logic to handle the partial barriers.

There are lots of ways to proceed. I modified the visualization to illustrate how indexing the sample times can be used to restrict the barrier penetration test and select those paths that survive the first barrier.

Clear[PartialBarrierTest];
PartialBarrierTest[path_, iSt_, iFn_, U_, L_] :=
Module[{subPath = path[[iSt ;; iFn]], subMax, subMin},
{subMax, subMin} = {Max[subPath], Min[subPath]};
(And[Max[subPath] < U, Min[subPath] > L])]

Module[{S0 = 100, r = 0.01, \[Sigma] = 0.25, L = 80, U = 120,
T = 1, \[CapitalDelta]t = 0.01, nTimes = 100, nPaths = 10,
seed = 123, drift, diff, times, paths, remainingPaths, PathToLine2},
SeedRandom[seed];
times = Range[0, T, \[CapitalDelta]t];
PathToLine2[path_, iSt_, iFn_] :=
Line[Transpose[{times[[iSt ;; iFn]], path[[iSt ;; iFn]]}]]; {drift,
diff} = {Exp[(r - \[Sigma]^2/2) \[CapitalDelta]t],
diff = \[Sigma] Sqrt[\[CapitalDelta]t]};
paths = Table[GBMPathCompiled[S0, drift, diff, nTimes], {nPaths}];

With[{iSt = 21, iFn = 41},
remainingPaths =
Select[paths, PartialBarrierTest[#, iSt, iFn, U, L] &]];
Graphics[{{Black, PathToLine2[#, 1, 41] & /@ paths}, {Blue,
PathToLine2[#, 41, 101] & /@ remainingPaths}, {Red,
Line[{{0.2, L}, {0.4, L}}], Line[{{0.8, L}, {1, L}}],
Line[{{0.2, U}, {0.4, U}}], Line[{{0.8, U}, {1, U}}]}},
Axes -> Automatic, AspectRatio -> 0.25]]

• Hi Tom, I edited your post to format your code. You can do this by selecting the relevant text and clicking the {} button above the editing area. Also, this is really an edit to the existing post. Could you please move the content to your other answer and delete this one? Mar 7 '13 at 4:12
• Hello Verbia, I'm sorry but I don't understand how to comply with your request. I don't see a {} button. I'm not sure how to move content and/or delete an entry. So many sites on the web, so many rules and conventions. I just tried to answer a question. Mar 8 '13 at 6:14
• @Tom Gladd Thank you. I found a way to do implement it based on your snippet. Could you maybe give me a hint how to implement a function which calculates a probability and takes as input two consecutive values in your framework (in PartialBarrierTest`) I want to calculate the probability that the GBM leaves the barriers between two discrete time stamps (similar to BBridge) and as input the functions expects the two values. THX Mar 15 '13 at 8:07