def METHOD_NAME(n, w, alpha, beta):
    np.random.seed(1)
    y = np.zeros(n)
    sigma2 = np.zeros(n)
    p = len(alpha)
    q = len(beta)
    w_vals = w < 0
    alpha_vals = np.any(alpha < 0)
    beta_vals = np.any(beta < 0)
    if np.any(np.array([w_vals, alpha_vals, beta_vals])):
        raise ValueError("Coefficients must be nonnegative")
    if np.sum(alpha) + np.sum(beta) >= 1:
        raise ValueError(
            "Sum of coefficients of lagged versions of the series and lagged versions of volatility must be less than 1"
        )
    # initialization
    if q != 0:
        sigma2[0:q] = 1
    for k in range(p):
        y[k] = np.random.normal(loc=0, scale=1)
    for k in range(max(p, q), n):
        psum = np.flip(alpha) * (y[k - p : k] ** 2)
        psum = np.nansum(psum)
        if q != 0:
            qsum = np.flip(beta) * (sigma2[k - q : k])
            qsum = np.nansum(qsum)
            sigma2[k] = w + psum + qsum
        else:
            sigma2[k] = w + psum
        y[k] = np.random.normal(loc=0, scale=np.sqrt(sigma2[k]))
    return y