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.
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.