The formula that you gave is not correct, by just looking at the maximum lag. If you are really into finding the truth about SARIMA formulas, I recommend that you check out the book "Time Series Analysis: Forecasting and Control", by Box, Jenkins, and Reinsel. ARIMA models are also known as Box-Jenkins models, that is why.
Now talk about (2,1,1)(1,1,1)12 model. First let me get the name right to facilitate the explanation: (p=2,d=1,q=1)(P=1,D=1,Q=1)S=12.
There are totally six polynomials: two for AR part (involving p and P), two for MA part (involving q and Q), two for differencing part (involving d and D).
Start with the differencing part, which is the product of two polynomials: (1-B)(1-B^12). This means the series is preprocess by differencing: Z[t] = Y[t] - Y[t-1] - Y[t-12] + Y[t-13].
Then the series Z[t] is modeled by an SARMA model, notice the "I" in the middle is gone.
Now look at the AR part, which is the product of two polynomials: (1 - phi1* B - phi2* B^2)(1 - phi3* B^12). Expand it I get: (1 - phi1* B - phi2* B^2 - phi3* B^12 + phi1*phi3* B^13 + phi2*phi3*B^14).
Now look at the MA part, which is the product of two polynomials: (1-theta1 * B)(1- theta2*B^12). Expand it I get: (1-theta1*B-theta2*B^12+theta1*theta2*B^13).
Now put AR polynomial and MA polynomial in the equation:
(1 - phi1* B - phi2* B^2)(1 - phi3* B^12). Expand it I get: (1 - phi1* B - phi2* B^2 - phi3* B^12 + phi1*phi3* B^13 + phi2*phi3*B^14)Z[t] = (1-theta1*B-theta2*B^12+theta1*theta2*B^13)a[t].
Now put Z[t] on the left side of the equation, and all other terms on the right:
Z[t] = (phi1* B + phi2* B^2 + phi3* B^12 - phi1*phi3* B^13 - phi2*phi3*B^14)Z[t]+ (1-theta1*B-theta2*B^12+theta1*theta2*B^13)a[t]
Now apply back operator rules, and get:
Z[t] = (phi1* Z[t-1] + phi2* Z[t-2] + phi3* Z[t-3] - phi1*phi3* Z[t-13] - phi2*phi3*Z[t-14]+ a[t]-theta1*a[t-1]-theta2*a[t-12]+theta1*theta2*a[t-13]
Notice a[t] = 0. Z is the differenced series, which you should be able to calculate from Y series using the formula that I mentioned about. the remaining a values are innovations from the fitted model. So, you should be able to get Z[t] now.
Now look at Z[t] = Y[t] - Y[t-1] - Y[t-12] + Y[t-13] again. Re-arrange it, and I get Y[t] = Z[t] + Y[t-1] + Y[t-2] - Y[t-13]. So I can get Y[t] given all previous Y's and Z[t].
The explicit formula of Y[t] is far more complicated than the formula that you provided.