Higher-order integrators are used, but usually the way they perform the calculation is not through directly calculating higher order derivatives, but essentially through multiple force calculations. The Wikipedia page on symplectic integrators gives some information on this, including 3rd order and 4th order examples. I've personally used the 4th order integrator in semiclassical simulations, where we were concerned about high accuracy to preserve quantum phase interference effects (although I'm not sure it was actually necessary).
In condensed phase systems, we're often simulating at a constant temperature, frequently with a stochastic thermostat, so the error in question should be how well we preserve the ensemble, not the energy itself. The "BAOAB" Langevin splitting has been shown to have a coincidental cancellation of errors that makes it a higher order integrator at the cost of a lower order scheme.
Non-symplectic predictor-corrector schemes used to be used more frequently, but it's generally agreed that symplectic methods are superior (thanks to some guaranteed properties coming from them being canonical transformations). I've used 4th-order Adam-Bashforth-Moulton predictor-correctors, as well as a 6th order Gear predictor corrector (modified to use adaptive time steps; needed when I hit a system where fixed time steps didn't work). The Gear predictor-corrector is presented in an appendix of [EDIT: the first edition of] Allen and Tildesley's "Computer Simulation of Liquids," which for a long time was one of the best books to learn MD (it's a little out of date now, being as it was published in 1987, [EDIT: and it looks like the new edition quite reasonably removed this appendix]). The Gear PC tries to actually estimate the values of the derivatives (whereas ABM uses a higher discretization in time). Some of my own old notes on ABM integrator and 4th order symplectic integrators can be found under the heading "Qual Prep" here (warning -- you'll have to decipher my handwriting).
Overall, although higher-order integrators have certainly been used, I think the overall sense in the community is that the cost-benefit ratio just isn't worth it.
Also, how would one actually go about calculating the third derivative of position with respect to time?
This should be feasible just by applying the chain rule. However, I think this would indeed require the whole Hessian matrix, which rules it out for large-scale MD.