We probe ultra-low-frequency gravitational waves ( GWs ) with statistics of spin-down rates of milli-second pulsars ( thereafter MSPs ) by a method proposed in our prevous work ( Yonemaru et al . 2016 ) . The considered frequency range is 10 ^ { -12 } { Hz } \lesssim f _ { GW } \lesssim 10 ^ { -10 } Hz . 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 } ~ { } { sec } ^ { -1 } and \dot { h } < 8.1 \times 10 ^ { -18 } ~ { } { sec } ^ { -1 } in the directions of the Galactic Center and M87 , respectively . Approximating the GW amplitude as \dot { h } \sim 2 \pi f _ { GW } h , the bounds translate into h < 3 \times 10 ^ { -8 } and h < 4 \times 10 ^ { -8 } , respectively , for f _ { GW } = 1 / ( 1000 ~ { } { yr } ) . Finally , we give the implications to possible super-massive black hole binaries at these sites .