Finding Option Probability Density Using Local Volatility from Dupire Model This question is different than https://quant.stackexchange.com/questions/31050/pricing-using-dupire-local-volatility-model
I am reading about the Dupire local volatility model. I have found ways to calculate the local volatility so for my question we can assume that it is known.
More specifically, it is considered piecewise constant between strikes and tenors so I have a local volatility surface that should be defined for all times and strikes.
From here I am wondering how I can think about solving the Dupire equation and recover risk neutral probability densities. I am not super familiar with stochastic differential equations so I am hoping that I can receive help reasoning through the problem.
My attempt
The Dupire equation takes the form
$dS_t = \mu_t S_t dt + \sigma(S_t,t) S_t dW_t$ where $S_t$ is the stock price at time $t$, $\mu_t$ is the drift term, $\sigma$ is the local volatility and $W_t$ is a Wiener process. Additionally, $S_t|_{t=0} = S_0$.
For simplicity, take $\mu_t=0$.
It is at this point that I have been getting confused about how to define the necessary constraints to solve the Ito integral.
First, I take it that $S_0$ is the spot price of the underlying of the derivative. If we solve this SDE, do we then find the how the spot price evolves over time? How does that help with option pricing? If $S_0$ is not the spot price, what instead is it? Is it the price of an option with a given strike at $t=0$ and we are then solving for the option prices at that strike across time?
If the latter is the case, do I solve this SDE for each forward price I want to query and simply have my volatility function be a function of time?
Once I have priced my options, I simply plan on taking the second derivative of price with respect to strike at each tenor to find the risk neutral probability density. Assume for this problem that tenors and strikes are dense enough that differentiation makes sense.
I found all of this information in Lecture 1: Stochastic Volatility and Local Volatility by Jim Gatheral, http://web.math.ku.dk/~rolf/teaching/ctff03/Gatheral.1.pdf
 A: In Dupire's local volatility model, the risk-neutral density $\psi (t,x)$ and the local volatility surface $\sigma (t,x)$ are related by the Fokker-Planck PDE: $$\frac{\partial}{\partial t} \psi (t,x) = \frac{1}{2} \frac{\partial^2}{\partial x^2} (\sigma^2(t,x) x^2 \psi (t,x)) \qquad (\star )$$
Or, equivalently, $$\frac{\partial}{\partial t} \psi = L_t^* \psi$$ where $L_t^*$ is the adjoint of the infinitesimal generator of $S_t$, where $$dS_t = \sigma (t, S_t) S_tdW_t$$
Note that we have assumed that interest rates are zero (i.e. $r = 0$).
Proof: For simplicity, suppose $\sigma$ is bounded and Lipschitz in the second variable; this implies that the SDE for $S_t$ admits a unique strong solution. Moreover, boundedness of $\sigma$ gives us that $S_t$ is a square-integrable martingale (via an argument using Grönwall's inequality and Doob's maximal inequality).
If, in addition, $\sigma$ is smooth and bounded away from zero, $S_t$ admits a risk-neutral density supported on a subset of $(0,\infty)$, say $\psi (t,x)$, so that for any $f \geq 0$: $$E(f(S_t)) = \int_0^\infty f(x)\psi (t,x) dx$$
To show $(\star )$ we now assume that $f \in C_c^\infty(0,\infty )$, so that Itô's formula gives us: $$f(S_t) = f(S_0) + \int_0^t f^\prime (S_s) dS_s + \frac{1}{2} \int_0^t f^{\prime \prime} (S_s) \sigma^2(s, S_s)S_s^2 ds$$
Taking expectations in two different ways and equating them,
$$
\begin{align*}
\int_0^\infty f(x)\psi (t,x) dx &= E(f(S_t)) \\
&= f(S_0) + \frac{1}{2} \int_0^t E[f^{\prime \prime} (S_s) \sigma^2(s, S_s)S_s^2] ds \\
&= f(S_0) + \frac{1}{2} \int_0^t \int_0^\infty f^{\prime \prime} (x) \sigma^2(s, x)x^2 \psi(s,x) dx ds \\
&= f(S_0) + \frac{1}{2} \int_0^t \int_0^\infty f (x) \frac{\partial^2}{\partial x^2} \left[\sigma^2(s, x)x^2 \psi(s,x) \right] dx ds
\end{align*}$$
where in the last step we have integrated by parts twice and used the fact that $f$ is compactly supported in $(0,\infty )$.
Differentiating with respect to $t$ and using Leibniz's rule,
$$\int_0^\infty f(x) \frac{\partial}{\partial t}  \psi (t,x) dx= \int_0^\infty f (x) \frac{1}{2} \frac{\partial^2}{\partial x^2} \left[\sigma^2(s, x)x^2 \psi(s,x) \right]dx$$
Since this equation holds for all $f \in C_c^\infty (0, \infty)$, $(\star)$ follows.
A: I misunderstood what it meant to solve an SDE. Each realization of $S_t$ that I find when I solve the SDE is a realization of the random variable $S_t$ and not the true terminal value of $S_t$. The way I find a PDF is by finding many instances of this random variable and finding their density at each time point.
I failed to grasp that solving the SDE a single time does not give the ONE AND ONLY value of $S_t$, rather that it is a possible value that occurs with some frequency. I can solve this using a monte carlo method and running it many times to get an idea of the pdf
