Here is a link to a paper I co-wrote which uses PyMC3 to solve linear models with residuals modelled by an autoregressive moving average or ARMA(p, q) distribution, that is, by any stationary distribution according to Wold’s Theorem.
This may be of interest to others as it allows values of p and q > 1, which do not seem to be implemented as native PyMC3 distributions at this point. Code is available at:
Note that these models require the derivation of reasonably informative prior distributions to converge (e.g. using statsmodels). This is likely because using p+q ARMA parameters is redundant - the number of parameters can be reduced to m=max(p, q+1) ( Chib, Siddhartha, and Edward Greenberg. “Bayes inference in regression models with ARMA (p, q) errors.” Journal of Econometrics 64.1-2 (1994): 183-206.)