  We propose an algorithm based on variational quantum imaginary time evolution
for solving the Feynman-Kac partial differential equation resulting from a
multidimensional system of stochastic differential equations. We utilize the
correspondence between the Feynman-Kac partial differential equation (PDE) and
the Wick-rotated Schr\"{o}dinger equation for this purpose. The results for a
$(2+1)$ dimensional Feynman-Kac system obtained through the variational quantum
algorithm are then compared against classical ODE solvers and Monte Carlo
simulation. We see a remarkable agreement between the classical methods and the
quantum variational method for an illustrative example on six and eight qubits.
In the non-trivial case of PDEs which are preserving probability distributions
-- rather than preserving the $\ell_2$-norm -- we introduce a proxy norm which
is efficient in keeping the solution approximately normalized throughout the
evolution. The algorithmic complexity and costs associated to this methodology,
in particular for the extraction of properties of the solution, are
investigated. Future research topics in the areas of quantitative finance and
other types of PDEs are also discussed.
