In this paper, we present a nonparametric estimation procedure for the multivariate Hawkes point process. The timeline is cut into bins and—for each component process—the number of points in each bin is counted. As a consequence of earlier results in Kirchner [Stoch. Process. Appl., 2016, 162, 2494–2525], the distribution of the resulting ‘bin-count sequences’ can be approximated by an integer-valued autoregressive model known as the (multivariate) INAR(p) model. We represent the INAR(p) model as a standard vector-valued linear autoregressive time series with white-noise innovations (VAR(p)). We establish consistency and asymptotic normality for conditional least-squares estimation of the VAR(p), respectively, the INAR(p) model. After appropriate scaling, these time-series estimates yield estimates for the underlying multivariate Hawkes process as well as corresponding variance estimates. The estimates depend on a bin-size and a support s. We discuss the impact and the choice of these parameters. All results are presented in such a way that computer implementation, e.g. in R, is straightforward. Simulation studies confirm the effectiveness of our estimation procedure. In the second part of the paper, we present a data example where the method is applied to bivariate event-streams in financial limit-order-book data. We fit a bivariate Hawkes model on the joint process of limit and market order arrivals. The analysis exhibits a remarkably asymmetric relation between the two component processes: incoming market orders excite the limit-order flow heavily whereas the market-order flow is hardly affected by incoming limit orders. For the estimated excitement functions, we observe power-law shapes, inhibitory effects for lags under 0.003 s, second periodicities and local maxima at 0.01, 0.1 and 0.5 s.