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.