ARIMA models in R are handled in a state space framework. See the help file for stats::arima
which includes the following section:
Fitting methods
The exact likelihood is computed via a state-space representation of the ARIMA process, and the innovations and their variance found by a Kalman filter. The initialization of the differenced ARMA process uses stationarity and is based on Gardner et al (1980). For a differenced process the non-stationary components are given a diffuse prior (controlled by kappa). Observations which are still controlled by the diffuse prior (determined by having a Kalman gain of at least 1e4) are excluded from the likelihood calculations. (This gives comparable results to arima0 in the absence of missing values, when the observations excluded are precisely those dropped by the differencing.)
Missing values are allowed, and are handled exactly in method "ML".
The Kalman filter allows for the likelihood to be computed exactly when missing values are present.
Model selection via forecast::auto.arima
does not use sample autocorrelations. It uses Akaike's Information Criterion which is based on the likelihood, so there is no problem in computing it when there are missing values present.