Individual Super-massive Black Hole Binaries in non-circular orbits are possible nanohertz gravitational wave (GW) sources for pulsar timing arrays. We develop an accurate and efficient approach to compute pulsar timing residuals due to GWs from inspiraling massive black hole binaries in relativistic eccentric orbits. Our approach employs a Keplerian-type parametric solution to model third post-Newtonian accurate precessing eccentric orbits while a novel semi-analytic prescription is provided to model the effects of quadrupolar order GW emission. These inputs lead to a semi-analytic prescription to model timing residuals from relativistic massive BH binaries inspiralling along arbitrary eccentric orbits. Additionally, we provide a f ully analytic approach to model timing residuals from BH binaries inspiraling along moderately eccentric orbits, influenced by Boetzel et al. 2017. Relevant Codes are being incorporated into Enterprise and TEMPO2 for searching the presence of such binaries in the PTA data-sets.