# Adding conditions to stochastic differential equations

Consider the following process

dt=0.001; s=1; tf=10;
f[x_, y_] := 2 - x^2;
g[x_, y_] := y - x y^2;
sol = RandomFunction[ItoProcess[{
\[DifferentialD]x[t] ==
f[x[t], y[t]] \[DifferentialD]t + \[DifferentialD]w[t],
\[DifferentialD]y[t] == g[x[t], y[t]] \[DifferentialD]t}, {x[t],
y[t]}, {{x, y}, {0, s}}, t,
w \[Distributed] WienerProcess[0, .1]], {0, tf, dt}];
ListLinePlot[sol, PlotRange -> All]


which produces the following result I would like to assign two conditions during the solution:

1. When x[t]==y[t] make y[t]=0

For this I've tried applying

\[DifferentialD]y[t] ==If[x[t]==y[t],0.5\[DifferentialD]t, g[x[t], y[t]]\[DifferentialD]t]

which did not work. It seems that this options is not valid using ItoProcess, so which one is?

1. Make an upper bound of 1.5 to x[t], such that in the proximity of the boundary, the noise values are bounded from above.

Two examples:

• Say that x[t]==1.5 and the noise can take values between [-0.1, 0.1], then I would like that in this specific case the noise can take values between [-0.1,0].

• If for example x[t]=1.495, I would like that the noise can take the values between [-0.1,0.05] etc.

Another way of treating this problem by forcing x[t]==1.5 whenever x[t]>1/5 (which brings us to the first part of the question). Regardless to this option, I would like the second part of the question to be addressed.

Two related questions:

Stochastic ODE Integration problems using RandomFunction

Boundary condition for stochastic differential equation

• Could you clarify condition 1? Do you want to set y[t]=0 or something with 0.5? May 5, 2019 at 19:45
• @ChrisK, set y[t]==0 May 5, 2019 at 20:05

You could use the same trick I used here -- manually adding the noise to NDSolve instead of using RandomFunction[ItoProcess].

σx = 0.1;
sol = NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]],
WhenEvent[Mod[t, dt] == 0,
{x[t] -> Min[1.5, RandomVariate[NormalDistribution[x[t], Sqrt[dt] σx]]]}],
WhenEvent[x[t] >= 1.5, {x[t] -> 1.5}],
x == 0, y == s}, {x, y}, {t, 0, tf}][];


The first WhenEvent periodically adds noise, but gives you the ability to keep x[t] < 1.5 with Min. The second WhenEvent is in case the deterministic dynamics makes x[t] > 1.5 between noise injections.

Checking the results:

Plot[Evaluate[{x[t], y[t], 1.5} /. sol, {t, 0, tf}] Note that there are two time steps in this approach: 1) a fixed time step dt between noise injections and 2) an automatic time step chosen by NDSolve. We can see this by ListLinePloting the InterpolatingFunction (a trick I learned from this answer by @MichaelE2).

ListLinePlot[x /. sol, PlotRange -> {{4.92, 4.93}, {1.48, 1.502}},
PlotStyle -> PointSize[0.01], PlotMarkers -> Automatic,
Epilog -> Line[{{4.92, 1.5}, {4.93, 1.5}}]] As you can see, NDSolve takes a few steps between noise injections (the vertical jumps). I don't think this is actually a problem, since any numerical solution to an SDE is an approximation that introduces an artificial fixed time step. I strongly suspect that if you make dt smaller, the solution will converge on a true trajectory, but this is not my area of expertise.

If the existence of two distinct time steps bothers you, then you can equate them by using a FixedStep method.

sol = NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]],
WhenEvent[Mod[t, dt] == 0,
{x[t] -> Min[1.5, RandomVariate[NormalDistribution[x[t], Sqrt[dt] σx]]]}],
WhenEvent[x[t] >= 1.5, {x[t] -> 1.5}],
x == 0, y == s}, {x, y}, {t, 0, tf},
StartingStepSize -> dt, Method -> {"FixedStep", Method -> "ExplicitEuler"}
][];

ListLinePlot[x /. sol, PlotRange -> {{2.53, 2.54}, {1.48, 1.502}},
PlotStyle -> PointSize[0.01], PlotMarkers -> Automatic,
Epilog -> Line[{{2.53, 1.5}, {2.54, 1.5}}]] You can see there are no intermediate steps between noise injections. This should be the Euler-Maruyama method.

I'm afraid I don't have time to figure out how to have both automatic step sizes and this flexible way to bound the solution. It'd be great if Wolfram would add WhenEvent to RandomFunction[ItoProcess].

If you want noise on both equations, just add an extra action to the WhenEvent.

I skipped over the first criterion, to set y[t] -> 0 when x[t] == y[t]. @Xminer gave one solution in their answer. Here's another that doesn't require a loose definition of equality:

sol = NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]],
WhenEvent[Mod[t, dt] == 0, {
xold[t] -> x[t],
x[t] -> Min[1.5, RandomVariate[NormalDistribution[x[t], Sqrt[dt] σx]]],
y[t] -> If[(xold[t] < y[t] && x[t] > y[t]) || (xold[t] > y[t] && x[t] < y[t]), 0, y[t]]
}],
WhenEvent[x[t] >= 1.5, {x[t] -> 1.5}],
WhenEvent[x[t] == y[t], {y[t] -> 0}],
x == 0, y == s, xold == 0}, {x, y}, {t, 0, tf},
DiscreteVariables -> {xold}][];

Plot[Evaluate[{x[t], y[t], 1.5} /. sol], {t, 0, tf}] This watches for three ways x[t] == y[t] is possible: they're equal between noise events or at a noise event x[t] changed between x[t] < y[t] and x[t] > y[t] or vice versa.

• Good luck with the RandomFunction. Could you be more specific about how you want the noise near the boundary? AFAIK, the noise in numerically solving an SDE is not drawn from a uniform distribution as you suggest in your post. May 5, 2019 at 11:22
• I think my implementation here with Min achieves a truncated normal distribution as you want. May 5, 2019 at 11:38
• @CATrevillian It separates the deterministic part from the stochastic part, which basically follows the Euler-Maruyama method where noise is treated as normally distributed with variance Sqrt[dt] σx. May 5, 2019 at 12:59
• I think this answer is the best. (Otherwise in Table + If + Euler method). Why do you stick to RandomFunction? when using arbitrary stochastic processes (e.g for complex option pricing), it's pretty tough to stick to a specific function. May 5, 2019 at 14:23
• @jarhead Yeah, it's hard to come up with a MWE that is simple but not too simple. I hope we see some other answers too. May 5, 2019 at 18:40

Dirty but works.

ClearAll["Global*"];
modifiedNoise[x_] :=
Block[{},
First@RandomVariate[
TruncatedDistribution[{-0.1, Max[0, 0.1 - Abs[1.5 - x]]},
NormalDistribution[0, .1*\[Sqrt]dt]], 1]];
run := Block[{}, dt = 0.02; s = 1; tf = 10; \[Epsilon] = .01;
f[x_, y_] := 2 - x^2;
g[x_, y_] := y - x y^2;

(*setting time index*)
numstep = IntegerPart[tf/dt];
timeindex = IntegerPart@Range[1, numstep];
Table[t[i] = 0 + dt*i, {i, timeindex}];
(*initial condition*)
t = 0;
x = 0;
y = s;

(*Explicit Forward Euler*)
Table[
newx =
x[i - 1] + dt*(f[x[i - 1], y[i - 1]]) + modifiedNoise[x[i - 1]];
newy = y[i - 1] + dt*(g[x[i - 1], y[i - 1]]);
If[Abs[newx - newy] > \[Epsilon],
x[i] = Min[1.5, newx];
y[i] = newy;,
x[i] = Min[1.5, newx];
y[i] = 0;]
, {i, timeindex}];

(*Result*)
Print[ListLinePlot[{Table[{t[i], x[i]}, {i, {0}~Join~timeindex}],
Table[{t[i], y[i]}, {i, {0}~Join~timeindex}]},
PlotLabels -> {"x[t]", "y[t]"}]]]; by the way,
when you add WhenEvent[Abs[x[t] - y[t]] <= 0.03, {y[t] -> 0}] to Chris's answer,
you got I didn't know this way, so her/his answer is better,I think.

• Nice! I wonder why y[t] tanks at t==1` compared to OP's and my simulations though. May 5, 2019 at 16:56
• the difference between us is about constraints of y(when x[t]==y[t],y[t]=0) and volatility of noise. From OP's definition of g[x,y],once y turned into 0 then it's clear that y keeps 0 or my mistake. May 5, 2019 at 17:02
• Ah I see what that first constraint means now. I might have to revisit my answer in light of that. May 5, 2019 at 17:25
• It was enough to add one line for your code. I did not know that NDSolve with WhenEvent could fully simulate SDE,really impressive! May 5, 2019 at 17:36
• thanks for that, I'll add it to my answer too May 5, 2019 at 19:03