  The paper focuses on pricing European-style options on several underlying
assets under the Black-Scholes model represented by a nonstationary partial
differential equation. The proposed method combines the Galerkin method with
$L^2$-orthogonal sparse grid spline wavelets and the Crank-Nicolson scheme with
Rannacher time-stepping. To this end, we construct an orthogonal cubic spline
wavelet basis on the interval satisfying homogeneous Dirichlet boundary
conditions and design a wavelet basis on the unit cube using the sparse tensor
product. The method brings the following advantages. First, the number of basis
functions is significantly smaller than for the full grid, which makes it
possible to overcome the so-called curse of dimensionality. Second, some
matrices involved in the computation are identity matrices, which significantly
simplifies and streamlines the algorithm, especially in higher dimensions.
Further, we prove that discretization matrices have uniformly bounded condition
numbers, even without preconditioning, and that the condition numbers do not
depend on the dimension of the problem. Due to the use of cubic spline
wavelets, the method is higher-order convergent. Numerical experiments are
presented for options on the geometric average.
