Thanks in advance for you interest,

I'm seeking advice if there are ways to speed up the computation of my simulation, without reducing accuracy.

I've been producing a simulation of an SDE of the Ornstein-Uhlenbeck process on a regular simplex. In particular, I look at a simplex with $m$ number of spatial dimensions, with the stochastic process $x$ starting at a uniform random location and evolving form there a tendency mean reversion. Further, the point of reversion changes halfway through the simulation.

Because of the unusual domain the best way to get at this is through simulating many replicates and averaging them together (in fact, I'm interested in a function $R(x)$ of the process).

However, I've had difficulty getting the standard error of the simulation low enough to for my purposes (finding min and max values over input parameter combinations). The levers I've been playing with are (1) number of replicates, (2) the time increment of simulation, (3) the SDE simulation method (e.g. Milstein).

What I'm seeking is advice on is if there are ways I can change my code to induce quicker computation, without sacrificing accuracy (e.g., reducing number of replicates makes things quicker but less accurate).

Thanks again for any input:

(** Computes \[CapitalEpsilon][R(X)] **)

(* \[Eta],\[Kappa],m,q,T are input parameters *)
(* \[Eta]: upper limit of a simplex side *)
(* \[Kappa]: rate of mean reversion *)
(* m: number of spatial dimensions *)
(* T: time limit *) 

(* \[CapitalDelta]t: the time increment of the SDE simulation *)
(* reps: number of replicate simulations *) 

\[ScriptCapitalE]R[\[Eta]_, \[Kappa]_, m_, q_, 
  T_, \[CapitalDelta]t_ : 0.01, reps_ : 10^4] :=
 Block[{xopt0, xopt1, F, R, x0, xs, xts, B, M, \[Sigma], \[Mu], proc, 
   RTs, t, \[Mu]0, \[Mu]1},

  (** produces random initial condition for each replicate **)
  x0 = Map[\[Eta]*Join[#, {1 - Total[#]}] &, 
    RandomVariate[DirichletDistribution[ConstantArray[1, {m}]], reps]];

  (** produces label of each spatial variable: X_1,X_2,...,X_m *)
  xs = Table[Symbol["x" <> ToString[i]], {i, m}];
  xts = Map[(#[t] &), xs];


  (** produces the deterministic dynamics and noise correlation \
matrices  **)
  (* mean reversion point is called xopt - 
  there are two because dynamics change halfway *)
  xopt0 = ConstantArray[\[Eta]/m, {m}]; 
  xopt1 = 
   Clip[If[q < 1.00, 
     Prepend[ConstantArray[\[Eta]/m - 
        1/m Log[(m - 1) q/(1 - q)], {m - 1}], \[Eta]/
       m + (m - 1)/m Log[(m - 1) q/(1 - q)]], 
     Prepend[ConstantArray[0, {m - 1}], \[Eta]]], {0, \[Eta]}]; 

  (* matrices constrain movement to simplex domain *)
  B[x_] := DiagonalMatrix[Map[Boole[# > 0] &, x]];
  M[x_] := Total[Map[Boole[# > 0] &, x]];

  (* \[Sigma] is the correlation matrix, \[Mu] the deterministic \
transition matrix *)
  \[Sigma] = 
   B[xts] . (IdentityMatrix[m] - 
      1/M[xts]*ConstantArray[1, {m, m}] . B[xts]);

  \[Mu][xoptin_] := \[Sigma] . (\[Kappa] (xoptin - xts));

  (* two deterministic matrices because dynamics change halfway *)
  \[Mu]0 = \[Mu][xopt0];
  \[Mu]1 = \[Mu][xopt1];

  (* provides function which random vector process X from SDE is fed \
into *)
  F[x_] := 1 - E^-x;
       R[xvec_] := 
   q*F[xvec[[1]]] + (1 - q)/(m - 1)*(Sum[F[xvec[[i]]], {i, 2, m}]);

  (** simulation of replicate trajectories **)
  (* halfway dynamics change *)
  (* 'proc' takes random initial condition and produces a trajectory *)
  (* trajectory is fed into 'R' *)
  proc[x0in_] :=
   Block[{X0, X1},
    X0 = 
     RandomFunction[
      ItoProcess[{\[Mu]0, \[Sigma]}, {xs, x0in}, {t, 0.}], {0., T/
       2, \[CapitalDelta]t}, Method -> "Milstein"];
    X1 = 
     RandomFunction[
      ItoProcess[{\[Mu]1, \[Sigma]}, {xs, X0["LastValues"][[1]]}, {t, 
        0.}], {T/2, T, \[CapitalDelta]t}, Method -> "Milstein"];
    TimeSeriesMap[R, TimeSeriesInsert[X0, X1]]["ValueList"]

    ];

  (* replicates are paralellized over random initial condition values \
*)RTs = Parallelize@Table[proc[x0[[i]]], {i, reps}];

  (* expected value over time and space is result *)
  Mean[Flatten[RTs]]

  ]

(* on the order of 3 minutes *)
AbsoluteTiming[\[ScriptCapitalE]R[1.00, .50, 3, 1.00, 1.00]]