For the first time, higher excited states of shallow acceptors up to the 8S and 5P states are calculated, using a method based on the Baldereschi-Lipari theory including the cubic correction. The eigenvalues and eigenvectors of the effective-mass Hamiltonian for shallow acceptor states were obtained by the finite-element method. The resulting sparse matrix is diagonalized by a newly developed method based on Arnoldis algorithm. Except for the lowest n, each hydrogenlike state nL splits into two levels when spherical spin-orbit coupling increases from 0 to 1. This results in crossing and repulsion of levels with different n. The spectra are thus shown to have totally different structure in the real acceptor regime 0.6 in contrast to exciton spectra for which 0.1. The calculated spectra are in agreement with available experimental data, especially in the case of higher excited states for which central-cell correction is negligible. The spectra of the shallow acceptors in ZnTe, CdTe, and InP are calculated and compared with the experimental ones.