This work is a continuation in a series of papers by the last two authors [cf., E. M. Stein and S. Wainger, Am. J. Math. 121, 1291–1336 (1999; Zbl 0945.42009); J. Anal. Math. 80, 335–355 (2000; Zbl 0972.42010)] pertaining to the study of discrete analogues of properties of singular integrals and maximal functions in harmonic analysis that depend, ultimately, on the geometry of the Euclidean space. The current paper considers a discrete analogue of the spherical average maximal operator \[ \sup_{\lambda>0} | f\ast d\sigma_\lambda| . \] This operator is bounded on \(L^p (\mathbb{R}^d)\) whenever \(p>d/(d-1)\) while \(d\geq 2\). Here \(d\sigma_\lambda\) denotes the normalized surface measure on the sphere of radius \(\lambda\). The discrete analogue of this operator takes the form \[ \sup_{\lambda>0}\frac{1}{N(\lambda)}\sum_{| m| =\lambda} f(n-m) \] where \(m,n\) take values in \(\mathbb{Z}^d\). This operator is bounded on a range of \(L^p\) spaces but, in contrast with the Euclidean case, it is necessary (and sufficient) that \(p>d/(d-2)\) while \(d\geq 5\) in order that this operator is bounded on \(\ell^p(\mathbb{Z}^d)\). As in previous results of the authors in this series of papers, considerations of Euclidean geometry like curvature need to be replaced by number theoretic concerns like class numbers. This is the principal reason for any differences in mapping properties. Analogues of ideas in sampling theory dating to Plancherel and Pólya also play a curious role here.