Stochastic volatility example with AR(1) instead of Gaussian random walk

I’ve been trying to implement and estimate, with pymc3, a basic stochastic volatility (SV) model of the following form:

r_t = exp{h_t/2}*e_t

h_t = r_0 + r_1*h_{t-1} + n_t

where r_t is the return process and h_t the (latent) log-variance process following a AR(1) process. My code (MWE) for this looks as follows:

import numpy as np
import pymc3 as pm

# simulate some random data
np.random.seed(13)
data = np.random.randn(10)

# SV model with AR
with pm.Model() as model:
    nu = 2
    rho = pm.Uniform("rho", -1, 1)

    h = pm.AR("h", rho=rho, sigma=1, shape=len(data))

    volatility_process = pm.Deterministic(
        "volatility_process", pm.math.exp(h / 2) ** 0.5
    )

    r = pm.StudentT("r", nu=nu, sigma=volatility_process, observed=data)
    prior = pm.sample_prior_predictive(10)
    # trace = pm.sample(10)

But running the above results in the following error message:

Traceback (most recent call last):
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py", line 801, in _draw_value
    return dist_tmp.random(point=point, size=size)
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\continuous.py", line 1979, in random
    point=point, size=size)
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py", line 638, in draw_values
    raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
ValueError: Cannot resolve inputs for ['Elemwise{mul,no_inplace}.0']

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 9, in <module>
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\sampling.py", line 1495, in sample_prior_predictive
    values = draw_values([model[name] for name in names], size=samples)
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py", line 620, in draw_values
    size=size)
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py", line 810, in _draw_value
    size=None))
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\continuous.py", line 1979, in random
    point=point, size=size)
  File "C:\Users\jrilla\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py", line 638, in draw_values
    raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
ValueError: Cannot resolve inputs for ['Elemwise{mul,no_inplace}.0']

Moreover, it is exactly the prior = ... line that causes the error. Note that I am using pm.AR() instead of pm.AR1(), but this would result in the same error. I don’t really understand why it does not work as expected. I am able to run the (simplified) SV example as provided in the pymc3 documentation:

# SV model with GaussianRandomWalk
with pm.Model() as model:
    nu = 2
    sigma = pm.Exponential("sigma", 1.0, testval=1.0)

    s = pm.GaussianRandomWalk("s", sigma=sigma, shape=len(data))

    volatility_process = pm.Deterministic(
        "volatility_process", pm.math.exp(-2 * s) ** 0.5
    )

    r = pm.StudentT("r", nu=nu, sigma=volatility_process, observed=data)
    prior = pm.sample_prior_predictive(10)
    # trace = pm.sample(10)

where they show the example for a Gaussian random walk (GRW), instead of the general AR process I want to use. Since the GRM is just a specific AR and it works for GRM, I don’t see why it shouldn’t also work for the general AR. As can be seen in the code, I basically just replace pm.GaussianRandomWalk(...) with pm.AR(...) (each with their required arguments). I am also able to simply implement/estimate the AR process itself:

# Simple AR
with pm.Model() as model:
    rho = pm.Uniform("rho", -1, 1)

    h = pm.AR("h", rho=rho, sigma=1, shape=len(data), observed=data)
    prior = pm.sample_prior_predictive(10)
    #trace = pm.sample(10)

which works fine as well, so I assume I am not making a mistake with defining the AR. It’s only when the AR is used as the latent process that the error arises.

Any idea on what the issue is here or what I’m doing wrong? I can’t seem to wrap my head around it.
Thanks!

Hi Jordi,
Thanks for reporting! It looks like a bug: a shape issue in sample_prior_predictive. If you comment this line out and uncomment trace = pm.sample(), your model should sample correctly.

@lucianopaz, are you aware of this problem? Is it worth filing a GitHub issue?
Here is the full traceback when trying to sample from the prior predictive:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    817                 try:
--> 818                     return dist_tmp.random(point=point, size=size)
    819                 except (ValueError, TypeError):

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
   2016         nu, mu, lam = draw_values([self.nu, self.mu, self.lam],
-> 2017                                   point=point, size=size)
   2018         return generate_samples(stats.t.rvs, nu, loc=mu, scale=lam**-0.5,

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    654             if to_eval == missing_inputs:
--> 655                 raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
    656             to_eval = set(missing_inputs)

ValueError: Cannot resolve inputs for ['Elemwise{mul,no_inplace}.0']

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-10-5e17a642d4eb> in <module>
     13 
     14     r = pm.StudentT("r", nu=nu, sigma=volatility_process, observed=data)
---> 15     prior = pm.sample_prior_predictive(10)
     16     trace = pm.sample(10)

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1848     names = get_default_varnames(vars_, include_transformed=False)
   1849     # draw_values fails with auto-transformed variables. transform them later!
-> 1850     values = draw_values([model[name] for name in names], size=samples)
   1851 
   1852     data = {k: v for k, v in zip(names, values)}

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    635                                         point=point,
    636                                         givens=temp_givens,
--> 637                                         size=size)
    638                     givens[next_.name] = (next_, value)
    639                     drawn[(next_, size)] = value

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    825                     with _DrawValuesContextBlocker():
    826                         val = np.atleast_1d(dist_tmp.random(point=point,
--> 827                                                             size=None))
    828                     # Sometimes point may change the size of val but not the
    829                     # distribution's shape

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
   2015         """
   2016         nu, mu, lam = draw_values([self.nu, self.mu, self.lam],
-> 2017                                   point=point, size=size)
   2018         return generate_samples(stats.t.rvs, nu, loc=mu, scale=lam**-0.5,
   2019                                 dist_shape=self.shape,

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    653         while to_eval or missing_inputs:
    654             if to_eval == missing_inputs:
--> 655                 raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
    656             to_eval = set(missing_inputs)
    657             missing_inputs = set()

ValueError: Cannot resolve inputs for ['Elemwise{mul,no_inplace}.0']

@jordyril, the problem is that neither the AR or the AR1 distributions have a random method. So you cannot sample from their prior using sample_prior_predictive. The only time series distribution that has a random method implemented is GaussianRandomWalk, so in practice, it’s the only one that you could use with sample_prior_predictive.

1 Like

Ooh, nice gotcha, thanks Luciano!
Is it possible to do prior predictive checks manually though (i.e without sample_prior_predictive)? Obviously not with PyMC random methods, as you explained, but like with numpy or scipy distributions?

I don’t think I understand what you mean regarding numpy and scipy. Are you asking if it’s possible to write a function that only used numpy and scipy, and generated samples from the prior of your model? Yes, you can do that but you’ll have to completely rewrite your model. Besides, while you rewrite it, you would have to actually implement our distributions missing random methods yourself. At that point, it would be more efficient to simply write the missing random method in pymc3.

Yes that’s exactly what I meant – I was trying to see if there was a quick workaround, but it doesn’t seem so. Plus, from a quick web search, I didn’t see AR distributions that you could sample from quickly in numpy and scipy (looks like you’d need statsmodels).
In any case, adding random methods to our timeseries distributions seems like a valuable addition!