Since the discovery of the first exoplanet around a Sun-like star in 1995, thousands of exoplanets have been detected in the past decades. However, none of them are like our Earth, which is a temperate rocky planet orbiting around a G star. Although a few high precision spectrographs were recently designed to detect sub-m/s radial velocity (RV) signals caused by Earth twins, the current RV modeling is not precise enough for the detection of such small signals. First, previous noise modeling does not choose optimal noise models through model comparison and thus suffer from overfitting or underfitting problems. Second, barycentric correction used in RV data reduction ignores the coupling between the Earth's motion and the stellar reflex motion, ignores the chromatic aberration and relativistic effects in the target system, and suffers from astrometric errors and errors in Earth's ephemeris. To avoid these problems, I have developed the Agatha and PEXO algorithms to mitigate correlated noise caused by stellar activity, instrument and atmosphere by choosing optimal noise models in the Bayesian framework, and to model the motions of the Earth and the target system as well as relativistic effects in them comprehensively. PEXO is able to model timing to nanosecond, astrometry to μas, and RV to μm/s. It will be used to synthesize astrometric, RV, timing and photometric data in order to detect Earth twins and exomoons. PEXO is also able to test general relativity in binary stars and to detect gravitational waves in pulsar timing data.