Supermassive black hole binaries, cosmic strings, relic gravitational waves from inflation, and first-order phase transitions in the early Universe are expected to contribute to a stochastic background of gravitational waves in the 10-9–10-7  Hz frequency band. Pulsar timing arrays (PTAs) exploit the high-precision timing of radio pulsars to detect signals at such frequencies. Here we present a time-domain implementation of the optimal cross-correlation statistic for stochastic background searches in PTA data. Due to the irregular sampling typical of PTA data as well as the use of a timing model to predict the times of arrival of radio pulses, time-domain methods are better-suited for gravitational-wave data analysis of such data. We present a derivation of the optimal cross-correlation statistic starting from the likelihood function, a method to produce simulated stochastic background signals, and a rigorous derivation of the scaling laws for the signal-to-noise ratio of the cross-correlation statistic in the two relevant PTA regimes: the weak-signal limit where instrumental noise dominates over the gravitational-wave signal at all frequencies, and a second regime where the gravitational-wave signal dominates at the lowest frequencies.