1
$\begingroup$

I'm trying to compute the average of the solution to an SDE by simulating some of its sample paths and then taking their Mean. The problem is my SDE is explosive and about 1 in 10 calls to RandomFunction will result in an Overflow message. Since I need ~10k paths for the LLN to take effect and give me a reliable average, I am stuck. Here is a minimal working example:

proc = ItoProcess[{{0}, {{x[t]^2}}, x[t]}, {{x}, {1}}, {t, 0}]
paths = RandomFunction[proc, {0., 5., 0.01}, 10];
list[t_] := paths[t][[2]][[2]];
me[t_] := Mean[list[t]];
Plot[me[t], {t, 0, 2}, AxesOrigin -> {0, 0}, AspectRatio -> 1, PlotRange -> {{0, 2} {0, 1}}]

Even with 10 paths it will often get stuck on some explosive path, and on 100 paths it almost certainly will. I would like a way of running the simulation and throwing away paths that result in an Overflow. I'm thinking of something along the lines of a For loop with a Throw-Catch in it, but couldn't find anything similar on the internet and am completely stuck.

I'm aware that this would result in slightly a biased sample; my preference would have been to get it to plot the solution stopped when its absolute value hits some threshold, but this seems unattainable since ItoProcess seems not to care what process I ask it to output: even if I substitute the first line above with

proc = ItoProcess[{{0}, {{x[t]^2}}, 1}, {{x}, {1}}, {t, 0}]

effectively asking it to output the constant 1, calls to RandomFunction will still crash with the same frequency. So I'll happily stick with just ignoring all the paths that result in an Overflow, unless of course there is a clever workaround.

Improve this question
$\endgroup$
4
  • $\begingroup$ "I'm aware that this would result in slightly a biased sample." Erm. I question that this would be only slightly biased if this happens once every ten times. Just image you would model a nuclear reaction in a civil nuclear power plant like this and what aftermath this could produce... $\endgroup$ – Henrik Schumacher Sep 30 '18 at 9:59
  • $\begingroup$ @HenrikSchumacher Fair enough. I didn't want to give the full details of my problem, but basically the reason why I'm not too concerned about this is that I'm only interested in the asymptotics of the expectation in zero. I actually proved that these are not affected by modifying the SDE outside of a neighbourhood containing the initial condition, so explosions don't really affect what I'm interested in. $\endgroup$ – Emilio Ferrucci Sep 30 '18 at 10:11
  • $\begingroup$ Then it would suffice to use a significantly smaller time horizon. For example, T = 0.1;paths = RandomFunction[proc, {0., T, T/50}, 10000]; works all fine for me. At least, this makes overflow less probable; maybe even impossible (reading from your comment, it seems that you have an argument that there is a positive minimal time until explosion). $\endgroup$ – Henrik Schumacher Sep 30 '18 at 10:15
  • $\begingroup$ @HenrikSchumacher Yes this is quite helpful, thank you. I guess I would have liked to have a way of systematically ignoring all explosions, and which allowed for an arbitrary time horizon. But this simple solution actually makes more sense theoretically. I can accept this as a solution if you post it. $\endgroup$ – Emilio Ferrucci Sep 30 '18 at 10:39
3
$\begingroup$

From you comment, I read off that you are actually interested only in the behavior close to $t = 0$. Then it would suffice to use a significantly smaller time horizon. For example, 

T = 0.1;
paths = RandomFunction[proc, {0., T, T/50}, 10000];

works all fine for me. At least, this makes overflow less probable; maybe even impossible (reading from your comment, it seems that you have an argument that there is a positive minimal time until explosion).

Improve this answer
$\endgroup$
2
$\begingroup$

I have not necessarily fully understood your question, but one approach is to specify a higher working precision for RandomFunction e.g.

proc = ItoProcess[{{0}, {{x[t]^2}}, x[t]}, {{x}, {1}}, {t, 0}]
paths = RandomFunction[proc, {0., 5., 0.01}, 10, WorkingPrecision -> 30];
list[t_] := paths[t][[2]][[2]];
me[t_] := Mean[list[t]];
Plot[me[t], {t, 0, 2}, AxesOrigin -> {0, 0}, AspectRatio -> 1, PlotRange -> {{0, 2} {0, 1}}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.