Waves observed in the photosphere and chromosphere of sunspots show complex dynamics and spatial patterns. The interpretation of high-resolution sunspot wave observations requires modeling of three-dimensional (3D) nonlinear wave propagation and mode transformation in the sunspot upper layers in realistic spot model atmospheres. Here, we present the first results of such modeling. We have developed a 3D nonlinear numerical code specially designed to calculate the response of magnetic structures in equilibrium to an arbitrary perturbation. The code solves the 3D nonlinear MHD equations for perturbations; it is stabilized by hyper-diffusivity terms and is fully parallelized. The robustness of the code is demonstrated by a number of standard tests. We analyze several simulations of a sunspot perturbed by pulses of different periods at a subphotospheric level, from short periods, introduced for academic purposes, to longer and realistic periods of 3 and 5 minutes. We present a detailed description of the 3D mode transformation in a non-trivial sunspot-like magnetic field configuration, including the conversion between fast and slow magneto-acoustic waves and the Alfvén wave, by calculation of the wave energy fluxes. Our main findings are as follows: (1) the conversion from acoustic to the Alfvén mode is only observed if the driving pulse is located out of the sunspot axis, but this conversion is energetically inefficient; (2) as a consequence of the cutoff effects and refraction of the fast magneto-acoustic mode, the energy of the evanescent waves with periods around 5 minutes remains almost completely below the level β = 1; (3) waves with frequencies above the cutoff propagate field aligned to the chromosphere and their power becomes dominating over that of evanescent 5 minute oscillations, in agreement with observations.