  The stochastic gravitational-wave background is imprinted on the times of
arrival of radio pulses from millisecond pulsars. Traditional pulsar timing
analyses fit a timing model to each pulsar and search the residuals of the fit
for a stationary time correlation. This method breaks down at
gravitational-wave frequencies below the inverse observation time of the array;
therefore, existing analyses restrict their searches to frequencies above 1
nHz. An effective method to overcome this challenge is to study the correlation
of secular drifts of parameters in the pulsar timing model itself. In this
paper, we show that timing model correlations are sensitive to sub-nanohertz
stochastic gravitational waves and perform a search using existing measurements
of pulsar spin-decelerations and pulsar binary orbital decay rates. We do not
observe a signal at our present sensitivity, constraining the stochastic
gravitational-wave relic energy density to $\Omega_\text{GW} ( f ) < 3.8 \times
10 ^{ - 9} $ at 450~pHz with sensitivity which scales as the frequency squared
until approximately 10 pHz. We place additional limits on the amplitude of a
power-law spectrum of $A_\star \lesssim 1.8\times10^{-14}$ for a reference
frequency of $f_* = 1~{\rm year} ^{-1} $ and the spectral index expected from
supermassive black hole binaries, $\gamma = 13/3$. If detection of a
supermassive black hole binary signal above 1 nHz is confirmed, this search
method will serve as a critical complementary probe of the dynamics of galaxy
evolution.
