Numerically Solving a 3d PDE with Stochastic Terms

I'm getting a bit confused if the procedure I'm doing is correct so any feedback would be great!

It's just a standard deterministic PDE for the price of a theoretic option, even if it's quite a messy one.

$S$ and $\sigma$ are stochastic and follow an SDE which I know. $t \in [0,T]$ and all other coefficients are constants.

I've done some research on solving 3d PDEs numerically and can't seem to find anything useful, let alone PDEs with stochastic components. So this is the procedure I've thought of:

1. Simulate paths of $S$ and $\sigma$ up to $T$. We need to inputting values of $S_0$ and $\sigma_0$ into the function to find the price of the option, and so we can't just work backwards by inputting terminal values of $S$ and $T$.
2. Evaluate $h(S_T, \sigma_T^2)$. This is the value of $C$ we have to start with.
3. Using finite difference terms, rearrange the PDE to be in the following form:
4. Numerical Scheme
5. For each iteration, n+1, of the scheme, set $\Delta_{S_{n+1}} = S_{n+1}-S_{n}$, $\Delta_{\sigma_{n+1}} = \sigma_{n+1}-\sigma_{n}$.
6. Compute $U(*1,*2,n)$ for $*1, *2=\{n-1,n,n+1\}$ and hence evaluate $U(m,k,n+1)$

Not sure if this is the right way to go about doing it so let me know if you have any comments/suggestions. Thanks in advance!

I personally like Finite Difference (FD) approaches, as they are quite suitable for this problems where the domain is generally rectangular, and you do not need too much work for evaluation of the Boundary Conditions.

What you wrote is a perfectly reasonable Explicit Euler (EE) scheme, and you can expect it to converge first order in time and second order in the "space variables", i.e. your scheme is a $O(\Delta t, \Delta S^2, \Delta\sigma^2)$.

Consider also that explicit schemes are easy to implement, but you lose stability. It might happen that if your $\Delta t$ is not "small enough" (see the Courant Number for the Heat Equation) your scheme will present spurious oscillations.

Knowing that your problem is linear, I would spend some time in investigating on some implicit schemes, e.g. Implicit Euler (IE), or the perfectly valid Alternating Direction Implicit (ADI) as @Canardini explained above.

In my experience, I've seen people implementing the Method Of Lines (MOL). That means that the time dimension is not discretised immediately. One firstly discretises the "space" variables, $S$ and $\sigma$, through you preferred technique (I recommend a 3 point technique such the one you proposed). Then you have an equation as:

$$\frac{dU(t, S_i, \sigma_j)}{dt} = f(t, U(t, S_i, \sigma_j), S_i, \sigma_j),$$

where $f$ is your spacial discretization. As you can see this is as easy as solving an ODE. If you're worried about accuracy then use $4^{th}$ order Runge-Kutta (RK4) or a Runge-Kutta-Felberg. If you're concerned about speed you might want to use a $2^{nd}$ order solver.

Consider that the EE and IE solver for this ODE coincides with the respective method for the original PDE.

Whatever you chose, I recommend to use some good numerical linear algebra library, as these kind of problems are described as "embarassingly parallelisable". Depending on the implementation language that you choose you have a multitude of options:

• MATLAB (integrated vectorization, but not free)

• Python: use numpy 1.11.1 as in here (precompiled with the Intel MKL libraries)

• Octave / Scilab (similar to MATLAB, but free)

If you're keener on programming: