Modeling MA(1) with heteroscedastic innovation term

Given a weekly time series, I am trying to model it with MA(1) process with heteroscedastic innovation term, with the variance hidden be some random walk

df = pd.read_csv('data/h1weekly.csv',
                 dtype={'IsCanceled': float})
df.index.freq = 'W-SUN'

series = pd.Series(np.log(df['IsCanceled'])).diff().dropna()
data = series.values

with pm.Model() as model:
    mu = pm.Normal('mu')
    theta = pm.Normal('theta')
    sigma = pm.GaussianRandomWalk('sigma', steps=len(data) - 1)
    epsilon = pm.Normal("epsilon", sigma=pm.math.abs(sigma))
    obs = pm.Normal('obs', mu=mu + theta * epsilon[:-1] + epsilon[1:], observed=data[1:])
    idata = pm.sample(1_000, chains=4)

The file is located here:

I would like to know:

In terms of your question about plotting, a great place to see how to do this can be seen in many of the example notebooks, PyMC Example Gallery — PyMC example gallery