Here’s a possible workaround…though I still prefer the type of solution involving a new multivariate class…
is it possible to define a Deterministic that’s a recursive function? this would be similar to the vol-update function in the GARCH class.
Consider the model again:
Y_t=\sqrt{V_{t-1}}\epsilon_t^s
log(V_t)=\alpha_v+\beta_vlog(V_{t-1})+\sigma_v\epsilon_t^v
\epsilon_t^s and \epsilon_t^v have correlation \rho
Set \epsilon^s=\rho\epsilon_t^v+\sqrt{1-\rho^2}\epsilon_t^n so that \epsilon_t^v and \epsilon_t^n are independent standard normals. Then we have
Y_t=\sqrt{V_{t-1}}[\rho\epsilon_t^v+\sqrt{1-\rho^2}\epsilon_t^n]
log(V_t)=\alpha_v+\beta_vlog(V_{t-1})+\sigma_v\epsilon_t^v
Now I can set up the model as follows
alpha=pm.Normal('alpha',0,sd=100)
beta=pm.Normal('beta',0,sd=10)
sigma_v=pm.InverseGamma('sigma_v',2.5,0.1)
epsilon_v=pm.Normal('error_v',mu=0,sd=1)
logV(t) is now a deterministic function of epsilon_v and the priors, but is recursive since it also depends on logV(t-1) this is the missing step
volatility=pm.Deterministic('volatility', pm.math.exp(.5*logV))
Yt=pm.Normal('Yt',mu=volatility*rho*error_v,sd=volatility*np.sqrt(1-rho*rho),observed=returns)
would something like this work?