For fun I'd like to use Hubble Astrometric data for example the observations listed at the bottom here to numerically estimate the orbit of 2014 MU69 as well as its uncertainties.
My plan is to use some combination of JPL Horizons, Hubble's TLEs, and Skyfield to get the J2000.0 position of the HST at the time of each exposure, and to get the position of the sun and major planets to generate the gravity field in which to integrate MU69's motion.
I understand I will have to retard the gravity from each source by its particular light-time, as well as correct for the light time for the HST images.
This will not be fast or efficient, it's strictly an exercise. At each time step I'll have to iterate and interpolate to find out, just for example "where would Jupiter have been in it's orbit such that it's gravity would be arriving right now".
I'd do that and calculate an initial orbit for MU69, then use that to calculate apparent positions for the HST data, calculate an error, then try a different starting state vector for MU69 and see if that's better or worse, and just use steepest descent to find a nominal orbit. From that, I can see how sensitive the fit is to various combinations of deviations from nominal.
I'm aware there may be smarter ways to do this, but to appreciate them it's better to do it brute force at least once. In the age of giga-flop laptops it's a viable option.
My question: Are there other things I need to consider?
Just for example, do I need to worry about time moving at different speeds at different distances from the sun (general relativity) or forces on MU69 besides gravity from the Sun and outer planets in order to get the level of accuracy relevant to comparison with the HST astrometry?
Again: ...in order to get the level of accuracy relevant to comparison with the HST astrometry — so I'm not looking for a list of arbitrarily small effects.