  GW170817 is a binary neutron star merger that exhibited a gravitational wave
(GW) and a gamma-ray burst, followed by an afterglow. In this work, we estimate
the Hubble constant ($H_0$) using broad-band afterglow emission and
relativistic jet motion from the Very Long Baseline Interferometry and Hubble
Space Telescope images of GW170817. Compared to previous attempts, we combine
these messengers with GW in a simultaneous Bayesian fit. We probe the $H_0$
measurement robustness depending on the data set used, the assumed jet model,
the possible presence of a late time flux excess. Using the sole GW leads to a
$20\%$ error ($77^{+21}_{-10}$ km/s/Mpc, medians, 16th-84th percentiles),
because of the degeneracy between viewing angle ($\theta_v$) and luminosity
distance ($d_L$). The latter is reduced by the inclusion in the fit of the
afterglow light curve, leading to $H_0=96^{+13}_{-10}$ km/s/Mpc, a large value,
caused by the fit preference for high viewing angles due to the possible
presence of a late-time excess in the afterglow flux. Accounting for the latter
by including a constant flux component at late times brings
$H_0=78.5^{+7.9}_{-6.4}$ km/s/Mpc. Adding the centroid motion in the analysis
efficiently breaks the $d_L-\theta_v$ degeneracy and overcome the late-time
deviations, giving $H_0 = 69.0^{+4.4}_{-4.3}$ km/s/Mpc (in agreement with
Planck and SH0ES measurements) and $\theta_v = 18.2^{+1.2}_{-1.5}$ deg. This is
valid regardless of the jet structure assumption. Our simulations show that for
next GW runs radio observations are expected to provide at most few other
similar events.
