What's the best way to implement an autoregressive moving average of a time series with a variable lag and window size?

Oh, set_subtensor is not an in-place operation, so your lagged function returns all nan values.

You cannot have missing values in your parameters. PyMC can do automatic interpolation if the data has missing values (i.e. c_data), but you have to define values for your model even in those cases, for example by assigning random variables to the missing values. The correct way to handle lagging is to create an (unobserved) initial state vector and prepend it to the data before convolution. See this notebook, cell 38 (“VAR as 2d convolution”) for an example.