I have an MA(4) process applied to the first order seasonal difference of $Y_t$ as follows:
$(1-B^s) Y_t = (1+\theta_1B+\theta_2B^2+\theta_3B^3+\theta_4B^4) Z_t$ where $Z_t \sim N(0,\sigma^2)$
This is equivalent to
$Y_t = \frac{(1+\theta_1B+\theta_2B^2+\theta_3B^3+\theta_4B^4)}{(1-B^s)} Z_t$ where $Z_t \sim N(0,\sigma^2)$
I understand how to refactor the denominator into an infinite order polynomial in the case when it has real roots - i.e. it can be decomposed into $(1-\phi B)$ terms which can then be inverted using the geometric sum identity $a/(1-r) = \sum_k ar^k$.
But in this case, because the denominator polynomial is a seasonal difference operator with S>2, the roots are complex (basically a unit circle on the complex plane).
So in short the question is how does one expand $(1-B^S)^{-1}$ into an infinite sequence? Or how do I determine the MA($\infty$) coefficients for the model above (numerically is fine, ie. happy to use numpy and truncate at some lag if it cannot be done analytically)