Recursive models are among the slowest in PyMC, especially when linear algebra operations are involved. It’s just the way it is for now. There might be a way to implement it in vector operations only; that’s why I wish we had started from there. I feel guilty for leading you down a needlessly general (and therefore needlessly difficult) path.
For example, looking at equation (2) from the paper:
\theta_t = \frac{\theta_{t-1}}{\gamma}\varepsilon_t
This is a Gaussian Random Walk in logs, i.e.:
\log \theta_t = \psi + \log \theta_{t-1} + \eta_t
Where \psi = -\log \gamma and \eta_t = \log \varepsilon_t.
So in PyMC, you could replace all the scanning with (all untested code):
gamma = pm.Beta('gamma')
psi = -pt.log(gamma)
#The .cumsum() makes it a GRW
log_theta = pm.Normal('log_theta', size=T).cumsum() + psi
The only hiccup is figuring out what to do with the innovations. The paper links the shape parameter \alpha to the predicted number of defaults in the previous period: \alpha_{t+1} = \gamma \alpha_t + N_{t-1} (hidden under equation 7). This would require a scan. This might be faster than what you’re doing now, though, because everything else is vectorized. Another option would be to abandon this dependency and just model \alpha in logs as a gaussian random walk. A third option would be to just abandon the time series component in the shape parameters all together, and just assume everything to do with time enters through \theta. I’d personally start with the 3rd option. In that case, you could just write:
epsilon = pm.Beta('epsilon', size=T)
theta = pm.Deterministic('theta', pt.exp(log_theta) * epsilon)
N = pm.Poisson('N', mu=theta, observed=data)
I’d probably start from there and see where it got me.