Find PDF of a stochastic process

Question

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.

Answer 1

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!