# Find PDF of a stochastic process

I consider the following Ito process

β = 2;
Ne = 100;
γ = 1;
process =
ItoProcess[{\[DifferentialD]s[
t] == -(β/Ne) s[t] i[t] \[DifferentialD]t -
Sqrt[β/Ne s[t] i[t]] \[DifferentialD]n1[
t], \[DifferentialD]i[
t] == (β/Ne s[t] - γ) i[t] \[DifferentialD]t +
Sqrt[β/Ne s[t] i[t]] \[DifferentialD]n1[t] -
Sqrt[γ i[t]] \[DifferentialD]n2[t]}, {s[t],
i[t]}, {{s, i}, {90, 10}},
t, {n1 \[Distributed] WienerProcess[],
n2 \[Distributed] WienerProcess[]}];


Is there a way to find the PDF of s and i or its discrete values for some discretization of t?

Update @Sjoerd Smit suggested to use the code

proc = ItoProcess[{1, 2}, {x, 0}, t];
rf = RandomFunction[proc, {0, 10}, 1000];
Plot[
PDF[SmoothKernelDistribution[rf["SliceData", 2]], x],
{x, -5, 10}
]


which works perfectly for some systems, but not on mine.

Then, @Josh Bishop pointed out that since I deal with a vector process, I should use

rf = RandomFunction[proc, {0., 5., 0.01}];
Plot3D[PDF[SmoothKernelDistribution[rf["SliceData", 2]], {s, i}], {s, 42, 49}, {i, 19, 27}]


However, the output is as follows, i.e., PDF is 0 everywhere.

• I updated my answer. Hopefully everything is clear :) – Sjoerd Smit Jan 29 at 10:31

For simple processes this can be done by using the function PDF:

proc = ItoProcess[{\[Mu], \[Sigma]}, {x, 0}, t];
PDF[proc[t], x]


E^(-((x - t [Mu])^2/(2 t [Sigma]^2)))/(Sqrt[2 [Pi]] Sqrt[ t [Sigma]^2])

If that doesn't work for your process, it's most likely that the analytic PDF simply cannot be computed. In that case you're pretty much down to sampling the process with RandomFunction to get a Monte Carlo approximation of the PDF. You can use a function like KernelMixtureDistribution or SmoothKernelDistribution to create a smooth PDF from the samples at a given time t.

For example, this is how to plot the approximate PDF at t = 2 for a simple process:

proc = ItoProcess[{1, 2}, {x, 0}, t];
rf = RandomFunction[proc, {0, 10}, 1000];
Plot[
PDF[SmoothKernelDistribution[rf["SliceData", 2]], x],
{x, -5, 10}
]


edit

Ah, I see the problem with this attempt. In the line

rf = RandomFunction[proc, {0., 5., 0.01}];


the process is sampled only once. You need multiple paths to approximate the PDF at a given point in time, so let's do 1000 samples (I'm going to assume that the process has already been defined as in the question):

rf = RandomFunction[process, {0, 5, 0.01}, 1000];


Let's say we want to inspect the PDF at t = 2.. Check the dimensions of the slice data:

t1 = 2.;
Dimensions[rf["SliceData", t1]]


{1000, 2}

First compute the domain of the PDF and pre-compute the distribution:

{{mins, maxs}, {mini, maxi}} = CoordinateBounds[rf["SliceData", t1]]
dist = SmoothKernelDistribution[rf["SliceData", t1]]


To plot the distribution, it is usually a good idea to keep the number of plot points in check by setting values for PlotPoints and MaxRecursion (the number of times Mathematica tries to refine the plot). My personal opinion is also that ContourPlot is almost always easier to interpret than Plot3D for these sort of things:

ContourPlot[
PDF[dist, {s, i}],
{s, mins, maxs},
{i, mini, maxi},
PlotPoints -> 20,
MaxRecursion -> 1,
PlotRange -> All
]


Hope this helps!

• Thank you very much, Sjoerd. I tried PDF, but it doesn't seem to work for my system or I do something wrong. If you have time, I would appreciate if you could demonstrate some uses of the functions you mention in your answer. Thank you! – Asatur Khurshudyan Jan 28 at 11:34
• @AsaturKhurshudyan If you made some new attempts based on my code, could you update the question with your new code? I hope I have time to give it some more attention later in the week, but right now I'm a bit strapped for time unfortunately. – Sjoerd Smit Jan 28 at 14:19
• @AsaturKhurshudyan You have a vector process so you'll need to use Plot3D. For example, you can try rf = RandomFunction[proc, {0., 5., 0.01}]; Plot3D[PDF[SmoothKernelDistribution[rf["SliceData", 2]], {s, i}], {s, 42, 49}, {i, 19, 27}]. – Josh Bishop Jan 28 at 19:03
• @AsaturKhurshudyan Yes, there is. You can use MarginalDistribution to convert the distribution into one of its maginals. E.g., PDF[MarginalDistribution[dist, 1], s] gives you the marginal over s and PDF[MarginalDistribution[dist, 2], i] the one over i. – Sjoerd Smit Jan 29 at 10:57
• Also, I noticed that the paths tend to start yielding complex values for larger values of t. You may want to think about what that means for the process, since it looks like something is ill-defined (unless s and i are in fact quantities that are allowed to be complex-valued, of course). – Sjoerd Smit Jan 29 at 11:00