How to solve the Black-Scholes equation using Feyman-Kac with two underlying assets? If we have the following Black-Sholes equation (PDE) for two independent stocks with the same volatility:
$$f_t +\sum^{d=2}_{i=1}rs_i\frac{\delta F}{\delta s_i} + \sigma^2\frac{1}{2}\sum^{d=2}_{i,j=1}\frac{\delta F}{\delta s_i,\delta s_j}-rf=0,$$
with terminal condition,
$$f(T,\boldsymbol{S}) = \max \left( \frac{1}{d}\sum^{d=2}_{i=1}s_i - K, 0 \right)$$
where $f$ is the option price, $s_i$ is price of stock $i$ at time $t$ and $\sigma$ is volatility, $K$ is exercise price of the optional at final time $T$, $r$ is the risk free interest rate, then what is the Feyman-Kac solution?
For example, in the corresponding one-dimensional case we have the solution:
$$f(t, s(t)) = \Bbb E \left[ e^{-r(T-t)} \max \left( s(t) - K, 0 \right) \right]$$
where $t$ is our starting time. In the two stocks case, do we just put $e^{-r(T-t)}$ in front of $\max(\frac{1}{d}\sum^{d=2}_{i=1}s_i -K, 0)$ as well? The reason I'm asking is that I want to use Monte Carlo simulations to estimate the price of the option in the two stocks case, or any $d>1$ case.
 A: Under the Black-Scholes framework, the answer is yes. However, the PDE that you provided does not exactly equal that from the B-S pricing model (missing $S^2$ terms in front of the second derivatives).
I will demonstrate this when the Brownian motions followed by the stock prices are uncorrelated (the argument is very similar otherwise). Suppose the SDEs under the risk-neutral martingale measure $\mathbb{Q}$ for the two assets are:
\begin{equation}
dX_t = X_t (rdt+\sigma dB_t),
\end{equation}
\begin{equation}
dY_t = Y_t (rdt+\sigma dW_t).
\end{equation}
Then, via a standard $\Delta$-hedging or self-financing argument, the PDE for $f(x,y,t)$ is:
\begin{equation}
f_t + r(xf_x+yf_y)+\frac{1}{2}\sigma^2(x^2f_{xx}+y^2f_{yy})-rf=0.
\end{equation}
Setting $g(x,y,t)=e^{-r(T-t)}f(x,y,t)$, the terminal condition remains the same whereas the new PDE is:
\begin{equation}
g_t + r(xg_x+yg_y)+\frac{1}{2}\sigma^2(x^2g_{xx}+y^2g_{yy})=0.
\end{equation}
So by Itô's lemma for the (hopefully) $C^2$ function $g(x,y,t) : \mathbb{R}^3 \to \mathbb{R}$, 
\begin{equation}
dg_t = \left( \frac{\partial g}{\partial t} + \frac{\sigma^2X_t^2}{2}\frac{\partial^2 g}{\partial x^2} + \frac{\sigma^2Y_t^2}{2}\frac{\partial^2 g}{\partial y^2} \right) dt + \dfrac{\partial g}{\partial x}dX_t + \dfrac{\partial g}{\partial y}dY_t = \\
\left( \frac{\partial g}{\partial t} + rX_t \frac{\partial g}{\partial x} + rY_t \frac{\partial g}{\partial y} + \frac{\sigma^2X_t^2}{2}\frac{\partial^2 g}{\partial x^2} +  \frac{\sigma^2Y_t^2}{2}\frac{\partial^2 g}{\partial y^2} \right) dt + \sigma X_t \dfrac{\partial g}{\partial x}dB_t + \sigma Y_t \dfrac{\partial g}{\partial y}dW_t.
\end{equation}
where all derivatives are evaluated at $(X_t,Y_t,t)$. The quadratic covariation term has vanished since two independent BMs are strongly orthogonal.
Now, integrate from $t$ to $T$ and take expectations, noting that the integrals with respect to the BMs are martingales and hence vanish in expectation. Moreover, the PDE for $g$ here is the same as the discounted Black Scholes PDE, hence the $dt$ term vanishes. We're therefore left with:
\begin{equation}
f(X,Y,t)=e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}_{t}[f(X_T,Y_T,T)|X_t=X,Y_t=Y].
\end{equation}
So plug in your terminal condition, simulate and discount - that's all. 
