  We probe ultra-low-frequency gravitational waves (GWs) with statistics of
spin-down rates of milli-second pulsars (MSPs) by a method proposed in our
prevous work (Yonemaru et al. 2016). The considered frequency range is
$10^{-12}{\rm Hz} \lesssim f_{\rm GW} \lesssim 10^{-10}$Hz, which cannot be
accessed by the conventional pulsar timing array. The effect of such
low-frequency GWs appears as a bias to spin-down rates which has a quadrupole
pattern in the sky. We use the skewness of the spin-down rate distribution and
the number of MSPs with negative spin-down rates to search for the bias induced
by GWs. Applying this method to 149 MSPs selected from the ATNF pulsar catalog,
we derive upper bounds on the time derivative of the GW amplitudes of $\dot{h}
< 6.2 \times 10^{-18}~{\rm sec}^{-1}$ and $\dot{h} < 8.1 \times 10^{-18}~{\rm
sec}^{-1}$ in the directions of the Galactic Center and M87, respectively.
Approximating the GW amplitude as $\dot{h} \sim 2 \pi f_{\rm GW} h$, the bounds
translate into $h < 3 \times 10^{-9}$ and $h < 4 \times 10^{-9}$, respectively,
for $f_{\rm GW} = 1/(100~{\rm yr})$. Finally, we give the implications to
possible super-massive black hole binaries at these sites.
