Abstract. The COS method for pricing European and Bermudan options with one underlying asset was developed in [F. Fang and C. W. Oosterlee, SIAM J. Sci. Comput., 31 (2008), pp. 826-848] and [F. Fang and C. W. Oosterlee, Numer. Math., 114 (2009), pp. 27-62]. In this paper, we extend the method to higher dimensions, with a multidimensional asset price process. The algorithm can be applied to, for example, pricing two-color rainbow options but also to pricing under the popular Heston stochastic volatility model. For smooth density functions, the resulting method converges exponentially in the number of terms in the Fourier cosine series summations; otherwise we achieve algebraic convergence. The use of an FFT algorithm, for asset prices modeled by Lévy processes, makes the algorithm highly efficient. We perform extensive numerical experiments. 1. Introduction. In financial markets traders deal in assets and options, like the well-known call and put options. Besides these, many "exotic" options have been defined that have more complex contract details and are not traded at regulated exchanges.One class of exotic option contracts is called the class of multicolor rainbow options, whose payoff may depend on multiple assets, like on the average or the maximum of asset prices. The value of the option depends on the contract details and on the underlying asset prices.Computational finance deals with numerical and computational questions regarding efficient option pricing and calibration. Usually, an asset price model is calibrated to liquidly available plain vanilla options (calls and puts) from a regulated exchange. For the valuation of the exotic options other computational methods are typically used. Option pricing techniques can be divided into the categories of Monte Carlo simulation, partial differential equation (PDE) methods, and Fourier-based methods. Often Monte Carlo methods are used to price high-dimensional option contracts. The method presented here can be seen as an alternative (deterministic) pricing technique, which can deal with multiasset option problems of medium-sized dimensionality, meaning two-dimensional (2D) to approximately five-dimensional (5D) integrals. The method we propose for pricing higher-dimensional options is based on the Fourier transform of the transitional density function and is especially suitable for asset price models in the class of Lévy processes.