Abstract This article considers fractionally integrated autoregressive moving-average time series models with conditional heteroscedasticity, which combines the popular generalized autoregressive conditional heteroscedastic (GARCH) and the fractional (ARMA) models. The fractional differencing parameter d can be greater than 1/2, thus incorporating the important unit root case. Some sufficient conditions for stationarity, ergodicity, and existence of higher-order moments are derived. An algorithm for approximate maximum likelihood (ML) estimation is presented. The asymptotic properties of ML estimators, which include consistency and asymptotic normality, are discussed. The large-sample distributions of the residual autocorrelations and the square-residual autocorrelations are obtained, and two portmanteau test statistics are established for checking model adequacy. In particular, non-stationary FARIMA(p, d, q)-GARCH(r, s) models are also considered. Some simulation results are reported. As an illustration,...