There’s no direct way to vectorize because it’s a recursive model, but you have options.
You could parametrize the simplex differently, such as using multiple gaussian random walk and then taking the softmax at each time step (you need one less random walk than the number of dimensions)
In that case the whole series can be vectorized easily because grw = normal(shape=t).cumsum()
Otherwise if you want to stick with Dirichlets you can use a Scan. That doesn’t vectorize but keeps the symbolic graph short: Time Series Models Derived From a Generative Graph — PyMC example gallery