We present a highly accurate, fully analytical model for the late inspiral, merger, and ringdown of black-hole binaries with arbitrary mass ratios and spin vectors and the associated gravitational radiation, including the contributions of harmonics beyond the fundamental mode. This model assumes only that nonlinear effects remain small throughout the entire coalescence, and is developed based on a physical understanding of the dynamics of late stage binary evolution, in particular on the tendency of the dynamical binary spacetime to behave like a linear perturbation of the stationary merger-remnant spacetime, even at times before the merger has occurred. We demonstrate that our model agrees with the most accurate numerical relativity results to within their own uncertainties throughout the merger-ringdown phase, and it does so for example cases spanning the full range of binary parameter space that is currently testable with numerical relativity. Furthermore, our model maintains accuracy back to the innermost stable circular orbit of the merger-remnant spacetime over much of the relevant parameter space, greatly decreasing the need to introduce phenomenological degrees of freedom to describe the late inspiral.Introduction.-Prior to the wide-ranging successes of numerical relativity (NR) that began with technical breakthroughs in 2005 [1-3] (see [4] for a recent review), the challenge of calculating the gravitational-wave emission from a pair of merging black holes was seen primarily as a problem on the boundary of nonlinear mathematics and computer science. The nonlinear nature of the partial differential equations describing general relativity was expected to manifest itself when the theory was pushed to describe the actual collision of black holes. The subsequent discovery that the radiation from the merger evolved very simply, smoothly connecting the amplitude and phase of the inspiral to those of the ringdown across all of the relevant parameter space, was a validation of two complementary efforts predicated on the assumed smallness of nonlinear effects throughout the entire coalescence -the close-limit approximation [5] culminating in the Lazarus project [6], and the Effective One-Body (EOB) approach [7,8]. However, although the smoothness of the merger has made it possible to create analytical models by extending post-Newtonian results to include free parameters, and tuning those to NR results [as is done in both EOB and the inspiral-merger-ringdown phenomenological (IMRPhenom) family of models [9]], there is currently no accurate model of the merger that is constructed analytically from first principles, rather than through a fit to NR.