Solution (source code)

= Solution

Differentiate the product of <Dirac delta functions> defining $\widetilde P$, regarding $B$ and $\sigma$ as independent density arguments. The random continuity equation is
$$
\partial_t\widetilde P=-\partial_B(\sigma B\widetilde P)+\frac1\tau\partial_\sigma(\sigma\widetilde P)-\partial_\sigma(f\widetilde P).
$$
For $0<s<t$, the causal <functional derivatives> follow from the explicit strain solution and $\widetilde B(t)=B_0\exp[\int_0^t\widetilde\sigma(u)\,du]$:
$$
\frac{\delta\widetilde\sigma(t)}{\delta f(s)}=e^{-(t-s)/\tau},\qquad \frac{\delta\widetilde B(t)}{\delta f(s)}=\widetilde B(t)\tau\left[1-e^{-(t-s)/\tau}\right].
$$
Consequently,
$$
\frac{\delta\widetilde P(t)}{\delta f(s)}=-e^{-(t-s)/\tau}\partial_\sigma\widetilde P-\tau\left[1-e^{-(t-s)/\tau}\right]\partial_B(\widetilde B(t)\widetilde P).
$$
In the <Furutsu–Novikov formula>, the covariance $\kappa\delta(t-s)$ samples the upper endpoint of the causal time integral. Symmetric regularization of the <white noise> gives half the mass there. The $B$ response vanishes as $s\uparrow t$, while the strain response tends to one, so $\mathbb E[f(t)\widetilde P(t)]=-(\kappa/2)\partial_\sigma P$. Averaging the continuity equation yields the <Fokker-Planck equation>
$$
\boxed{\partial_tP=-\partial_B(\sigma BP)+\frac1\tau\partial_\sigma(\sigma P)+\frac\kappa2\partial_\sigma^2P.}
$$
There is no direct diffusion term in $B$: its equation contains the time integral of <colored noise>, whereas the strain equation receives <Gaussian white noise>. The same result follows from the <Itô diffusion> with drift $(\sigma B,-\sigma/\tau)$ and diffusion vector $(0,\sqrt\kappa)$.