There are different methods to decide on the order of integration for a nonseasonal AR(I)MA model. Hyndman & Khandakar (2008, section 3.1) give pointers to the most commonly encountered ones. The most common type would be unit root tests, especially the Dickey-Fuller test, which Hyndman & Khandakar counsel against, since it biases towards more rather than fewer differences. Instead, they use a KPSS test (Kwiatkowski et al., 1992): you test for a unit root; if the test is significant, you difference and test again, until the test is not significant any more.
Yes, these are not the most recent papers, but
auto.arima() in the
forecast package for R still uses this approach, and that is pretty much as close to the gold standard in time series analysis as you can get.
After you have decided on the order of integration, you need to decide on AR and MA orders. Parsing ACF/PACF plots of successive residuals is the older Box-Jenkins approach; the more modern way would be to minimize an information criterion like the AICc. See the fuller description of how
auto.arima() decides on a model order and estimates.
In the present case,
auto.arima() would go for an ARIMA(1,1,1) model:
births <- read.table("daily-total-female-births-CA.csv",header=TRUE,sep=",",colClasses=c("Date","numeric"))
births_ts <- ts(births$births,frequency=365,start=births$date)
Since you work in Python, you may be interested in
pmdarima and in this SO thread: auto.arima() equivalent for python.