Time domain analysis of nonlinear systems with hysteretic damping is conducted. Specifically, the viscoelastic model proposed by Biot is examined. This hysteretic element represents an integral transform in the time domain. Thus, it yields integro-differential equations when it is incorporated into the system dynamics models. Two numerical methods are proposed to solve these equations. The first method approximates the kernel of this integral transform by a sum of exponentials making the computational cost minimal. The second method uses digital filters designed to match the transfer function, real and imaginary parts, of the Biot hysteretic element. These techniques are employed in calculating the response of a single-degree-of-freedom (SDOF) system with hysteretic damping and nonlinear stiffness subjected to deterministic, seismic, and random excitation. The method of statistical linearization is used to estimate the variance of the response of the SDOF system subjected to white noise. The accuracy of the results is verified by pertinent Monte Carlo studies. The presented approaches can be extended to treat multi-degree-of-freedom (MDOF) systems with hysteretic behavior.