This paper concerns the numerical solution of the two-dimensional time-dependent partial integro-differential equation that holds for the values of European-style options under the two-asset Kou jump-diffusion model. A main feature of this equation is the presence of a nonlocal double integral term. For its numerical evaluation, we extend a highly efficient algorithm derived by Toivanen in the case of the one-dimensional Kou integral. The acquired algorithm for the two-dimensional Kou integral has an optimal computational cost: the number of basic arithmetic operations is directly proportional to the number of spatial grid points in the semidiscretization. For effective discretization in time, we study seven contemporary implicit–explicit and alternating-direction implicit operator splitting schemes. All these schemes allow for a convenient, explicit treatment of the integral term. We analyze their (von Neumann) stability. Through ample numerical experiments for put-on-the-average option values, we investigate the actual convergence behavior as well as the relative performance of the seven operator splitting schemes. In addition, we consider the Greeks Delta and Gamma.