A search for a massive resonance decaying into a pair of standard model Higgs bosons, in a final state consisting of two b quark–antiquark pairs, is performed. A data sample of proton–proton collisions at a centre-of-mass energy of 13 TeV is used, collected by the CMS experiment at the CERN LHC in 2016, and corresponding to an integrated luminosity of 35.9. The Higgs bosons are highly Lorentz-boosted and are each reconstructed as a single large-area jet. The signal is characterized by a peak in the dijet invariant mass distribution, above a background from the standard model multijet production. The observations are consistent with the background expectations, and are interpreted as upper limits on the products of the s-channel production cross sections and branching fractions of narrow bulk gravitons and radions in warped extra-dimensional models. The limits range from 126 to 1.4 fb at 95% confidence level for resonances with masses between 750 and 3000 GeV, and are the most stringent to date, over the explored mass range.