The proximate problem is that the quantity a[t-1] + phi*b[t-1] + s[t-slen] is an aesara symbolic tensor, not a number, so numpy isn’t letting you store it in f. But the ultimate problem is that looping over the data this way isn’t going to give you what you want, because you aren’t taking into account the “innovations” at each time step, which all need to be estimated.
Copying from Forecasting: Principals and Practice, it looks like your model is an additive seasonal model with dampened trend, something like:
\begin{align} y_t &= \ell_{t-1} + \phi b_{t-1} + s_{t-m} + \epsilon_t \\ \ell_t &= \ell_{t-1} + \phi b_{t-1} + \alpha \epsilon_t \\ b_t &= \phi b_{t-1} + \beta \epsilon_t \\ s_t &= s_{t-m} + \gamma \epsilon_t \end{align}
Which can be put into state space form. For concreteness, let’s say the seasonal lag m=3, otherwise I have to type out a bunch of dots and stuff in the matrices.
A state space model is written:
\begin{align}x_{t+1} &= T x_t + R \epsilon_{t+1} \\ y_t &= Z x_t + H \eta_t \end{align}
with:
x_t = \begin{bmatrix} \epsilon_t \\ \ell_t \\ b_t \\ s_t \\ s_{t-1} \\ s_{t-2} \end{bmatrix}, T = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 &0 \\ \alpha & 1 & \phi & 0 & 0 & 1\\ \beta & 0 & \phi & 0 & 0 & 0 \\ \gamma & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1& 0 & 0 \end{bmatrix}, R = \begin{bmatrix} 1 \\ 0\\ 0\\ 0 \\ 0 \\ 0 \end{bmatrix}, Z = \begin{bmatrix} 1 & 1 & \phi & 1 & 0 & 0 \end{bmatrix}, H = 0.
This is a nice use-case for my PyMC Statespace “package” (ok it’s just a repo).
Statespace models estimated with a Kalman filter are inherently Bayesian; you always get confidence intervals thanks to the multivariate normal being conjugate with itself. You could also estimate this using statsmodels.api.tsa.statespace.ExponentialSmoothing. Your errors would be underestimated because they won’t account for model uncertainty (i.e., uncertainty about the parameter values themselves), but it’s still quite nice. My thing isn’t really production ready.