By the way, if you are interested in timeseries in PyMC, we have access to a nifty new library Aeppl that among other things, tries to infer the logp of arbitrary Scans that represent the generative process of a timeseries (i.e., the random draws / the graph implemented in rv_op): GitHub - aesara-devs/aeppl: Tools for an Aesara-based PPL.
I tried to use that when refactoring the AR distribution and it did work. We ended up using the manual vectorized implementation because the generated Scan was too slow for now (something that the folks at Aesara are trying to iron out). Anyway you might like to try it out. I have a notebook showing it within a PyMC model (warning: it has some black magic to glue things together that may no longer be needed / working):