Efficient GeometricBrownianMotionProcess Monte Carlo simulation

Question

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[#, #[[1]] >= Ti[[i]] && #[[1]] <= (Ti[[i]] + per) &] & /@ 
        {data["Paths"][[j]]}, {i, 1, Length[Ti]}];
      data2 = Flatten[data1, 2];
      value += 
        If[Max[Table[data2[[i]][[2]], {i, 1, Length[data2]}]] <= U && 
           Min[Table[data2[[i]][[2]], {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=

Answer 1

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

Roll your own GBM

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.

Answer 2

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]]