Fractional ARIMA simulation

Question

I'm trying to design from scratch a function to simulate a fractional ARIMA time series. Unfortunately I dont know how to add the fractional part to the function, while the AR and MA part are already implemented.

At the moment what I did is the ARMA, how do integrate the fractional differentiation part?

My code is:

import numpy as np

def ARMA(phi, theta, n):
    """
    Generate ARMA(p, q) time series data of length n.

    Parameters:
        phi (list): Coefficients of autoregressive terms.
        theta (list): Coefficients of moving average terms.
        n (int): Length of the time series.

    Returns:
        np.array: Generated ARMA(p, q) time series data.
    """
    p = len(phi)
    q = len(theta)

    # Generate white noise
    white_noise = np.random.normal(size=n)

    # Initialize the time series data with zeros
    arma_series = np.zeros(n)

    # Generate ARMA(p, q) time series
    for i in range(max(p, q), n):
        ar_term = np.dot(phi, arma_series[i-p:i][::-1])
        ma_term = np.dot(theta, white_noise[i-q:i][::-1])
        arma_series[i] = ar_term + ma_term + white_noise[i]

    return arma_series

# Example usage:
# Define ARMA coefficients
phi = [0.5, -0.2]
theta = [0.3, 0.4]
n = 100

# Generate ARMA time series data
arma_data = ARMA(phi, theta, n)
print(arma_data)