Exaggerated traits involved in species interactions have long captivated the imagination of evolutionary biologists and inspired the durable metaphor of the coevolutionary arms race. Despite decades of research, however, we have only a handful of examples where reciprocal coevolutionary change has been rigorously established as the cause of trait exaggeration. Support for a coevolutionary mechanism remains elusive because we lack generally applicable tools for quantifying the intensity of coevolutionary selection. Here we develop an approximate Bayesian computation (ABC) approach for estimating the intensity of coevolutionary selection using population mean phenotypes of traits mediating interspecific interactions. Our approach relaxes important assumptions of a previous maximum likelihood approach by allowing gene flow among populations, variable abiotic environments, and strong coevolutionary selection. Using simulated data, we show that our ABC method accurately infers the strength of coevolutionary selection if reliable estimates are available for key background parameters and ten or more populations are sampled. Applying our approach to the putative arms race between the plant Camellia japonica and its seed predatory weevil, Curculio camelliae, provides support for a coevolutionary hypothesis but fails to preclude the possibility of unilateral evolution. Comparing independently estimated selection gradients acting on Camellia pericarp thickness with values simulated by our model reveals a correlation between predicted and observed selection gradients of 0.941. The strong agreement between predicted and observed selection gradients validates our method.