  We present a python based parameter inference system for the gravitational
wave (GW) measured in the millihertz band. This system includes the following
features: the GW waveform originated from the massive black hole binaries
(MBHB), the stationary instrumental gaussian noise, the higher-order harmonic
modes, the full response function from the time delay interferometry (TDI) and
the gaussian likelihood function with the dynamic nested parameter sampler. In
particular, we highlight the role of higher-order modes. By including these
modes, the luminosity distance estimation precision can be improved roughly by
a factor of 50, compared with the case with only the leading order
($\ell=2,|m|=2$) mode. This is due to the response function of different
harmonic modes on the inclination angle are different. Hence, it can help to
break the distance-inclination degeneracy. Furthermore, we show the robustness
of testing general relativity (GR) by using the higher-order harmonics. Our
results show that the GW from MBHB can simultaneously constrain four of the
higher harmonic amplitudes (deviation from GR) with a precision of
$c_{21}=0.54^{+0.61}_{-0.82}$, $c_{32}=-0.65^{+0.22}_{-0.08}$,
$c_{33}=0.56^{+0.60}_{-0.76}$ and $c_{44}=1.57^{+2.34}_{-1.90}$, respectively.
