The model is defined by values for the AR and MA parameters (\(\phi\) and \(\theta\), respectively), along with the fractional differencing parameter d. When \(d\geq 0.5\), then the integer part is taken as \(m=\lfloor d+0.5\rfloor\), and the remainder (between -0.5 and 0.5) stored as d. For \(m=0\), the model is:
$$\left(1 - \sum_{i=1}^p \phi_i B^i\right)\left(1 - B\right)^d (y_t - \mu)=\left(1 + \sum_{i=1}^q \theta_i B^i\right) \epsilon_t$$
where B is the backshift operator (\(B y_t = y_{t-1}\)) and \(\epsilon_t\) is the innovation series. When \(m > 0\), the model is defined by:
$$y_t = (1 - B)^{-m}x_t$$
$$\left(1 - \sum_{i=1}^p \phi_i B^i\right)(1 - B)^d (x_t - \mu)=\left(1 + \sum_{i=1}^q \theta_i B^i\right) \epsilon_t$$
When stat.int = FALSE, the differencing filter applied to the innovations is not split into parts, and the series model follows the first equation regardless of the value of d. This means that \(\mu\) is added to the series after filtering and before any truncation. When stat.int = TRUE, \(x_t - \mu\) is generated from filtered residuals, \(\mu\) is added, and the result is cumulatively summed m times. For non-zero mean and \(m>0\), this will yield a polynomial trend in the resulting data.
Note that the burn-in length may affect the distribution of the sample mean, variance, and autocovariance. Consider this when generating ensembles of simulated data