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 fb$^{-1}$. 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.