We consider numerical methods for the efficient and stable solution of advanced, multidimensional partial integro-differential equations (PIDEs) and partial integro-differential complementarity problems (PIDCPs) arising in financial energy option valuation. Here the underlying uncertain factors, e.g. the electricity price, are modelled by exponential Lévy processes to account for the jumps that are often observed in the markets. These jumps give rise to the integral term in the PIDEs and PIDCPs and is nonlocal. For the effective numerical solution, we investigate in this project operator splitting methods. This broad class of methods has already been successfully applied and analysed in the special case of partial differential equations (PDEs). Many questions about their adaptation to PIDEs and PIDCPs are, however, still largely open. In this project we consider two main research topics: (1) infinite activity jumps and (2) swing options. The first topic concerns exponential Lévy processes with an infinite number of jumps in every time interval. The second topic deals with a popular type of exotic energy option that has multiple exercise times. We develop novel, second-order operator splitting methods of the implicit-explicit (IMEX) and alternating direction implicit (ADI) kind. We analyse their fundamental properties of stability, consistency, monotonicity and convergence. The acquired theoretical results are validated by ample numerical experiments.