Add GARCH11 model

[olafSmits]

Here's the GARCH, I hope @springcoil
This is my first pull request, so I hope I'm doing this right.

[springcoil]

Looks good to me. Now I understand what I was doing wrong. Will review
later in more depth.
On 20 Feb 2016 8:43 PM, "Olaf Smits" notifications@github.com wrote:

Here's the GARCH, I hope @springcoil https://github.com/springcoil

This is my first pull request, so I hope I'm doing this right.

You can view, comment on, or merge this pull request online at:

#987
Commit Summary

• Add GARCH11 model

File Changes

• M pymc3/distributions/timeseries.py
  https://github.com/pymc-devs/pymc3/pull/987/files#diff-0 (52)

Patch Links:

• https://github.com/pymc-devs/pymc3/pull/987.patch
• https://github.com/pymc-devs/pymc3/pull/987.diff

—
Reply to this email directly or view it on GitHub
#987.

[fonnesbeck]

Why is scan needed? I the function not vectorized/vectorizable?

[springcoil]

Basically it is not vectorizable, in other work (including my attempts this became apparent.) The original version by Olaf - I believe ported from R used a for loop, and this has been translated into a scan. Since scan is effectively the theano version of doing a for loop.

[olafSmits]

No, you can't vectorize it. The "definition" of vol[n+1] depends on vol[n], which in turn depends on vol[n-1], etc, all the way to the initial volatility.

[springcoil]

A problem is when I run my example I have an error produced.

J = 8
r = shared(np.array([28, 8, -3, 7, -1, 1, 18, 12]))
sigma1 = shared(np.array([15, 10, 16, 11, 9, 11, 10, 18]))
alpha0 = shared(np.array([10, 10, 16, 8, 9, 11, 12, 18]))

with Model() as garch:
    alpha1 = Normal('alpha1', 0, 1, shape=J)
    BoundedNormal = Bound(Normal, upper=(1 - alpha1))
    beta1 = BoundedNormal('beta1', 0, sd=1e6)
    mu = Normal('mu', 0, sd=1e6)


    obs = GARCH11('garchy_garch', initial_vol=alpha1, observed=r)

    step = NUTS()

My error is the following.
Traceback (most recent call last):
  File "pymc3/examples/garch_example.py", line 52, in <module>
    obs = GARCH11('garchy_garch', initial_vol=alpha1, observed=r)
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/pymc3-3.0-py3.5.egg/pymc3/distributions/distribution.py", line 25, in __new__
    return model.Var(name, dist, data)
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/pymc3-3.0-py3.5.egg/pymc3/model.py", line 256, in Var
    var = ObservedRV(name=name, data=data, distribution=dist, model=self)
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/pymc3-3.0-py3.5.egg/pymc3/model.py", line 520, in __init__
    self.logp_elemwiset = distribution.logp(data)
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/pymc3-3.0-py3.5.egg/pymc3/distributions/timeseries.py", line 118, in logp
    vol = self._get_volatility(x[:-1])
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/pymc3-3.0-py3.5.egg/pymc3/distributions/timeseries.py", line 114, in _get_volatility
    self.beta_1])
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/Theano-0.8.0.dev0-py3.5.egg/theano/scan_module/scan.py", line 365, in scan
    non_seqs.append(tensor.as_tensor_variable(elem))
  File "/Users/peadarcoyle/anaconda/envs/pymc3_examples/lib/python3.5/site-packages/Theano-0.8.0.dev0-py3.5.egg/theano/tensor/basic.py", line 208, in as_tensor_variable
    raise AsTensorError("Cannot convert %s to TensorType" % str_x, type(x))
theano.tensor.var.AsTensorError: ('Cannot convert None to TensorType', <class 'NoneType'>)

[olafSmits]

There's a couple of things going wrong with this example. The most important one is that I shouldn't have set the default values of the different parameters to None -- that's what's generating the errors. So in your example you're not passing sigma0 etc to GARCH11, and so it's just using the default None which throws your error.

I guess I should get rid of the default initialization to None. Should I just use some default distributions for the parameters when None is passed?

When I fix that I get errors thrown by scan. On my laptop scan complains if you use anything else then float32 due to casting from int64 or float64 to float32, so you need to make sure that all the input is in this dtype.

Another error is that there should be a lower bound on beta1 .
Finally, the parameters are all scalars. Only the returns r should be a vector.

Unfortunately, I still get TypeErrors that are associated with the NUTS() and I'm not sure how to fix it. It's something with the conversion between float32 and float64.

That said, I modified the example to this:

r = shared(np.array([28, 8, -3, 7, -1, 1, 18, 12, 15], dtype=np.float32))
sigma1 = shared(np.array(.2, dtype=np.float32))
alpha0 = shared(np.array(.5, dtype=np.float32))

with Model() as garch:
    alpha1 = Normal('alpha1', 0., 1., dtype='float32')
    BoundedNormal = Bound(Normal, lower=0., upper=(1 - alpha1))
    beta1 = BoundedNormal('beta1', 0.,sd=1e3, dtype='float32')

    obs = GARCH11('garchy_garch', omega=alpha0, alpha_1=alpha1,
                                beta_1=beta1, initial_vol=sigma1, observed=r,
                                dtype='float32')
    start = find_MAP(fmin=optimize.fmin_bfgs)

This code runs for me. Unfortunately it doesn't run with the NUTS().

[springcoil]

What do you run the trace on then - I'm trying it with Metropolis and getting errors.

I'm using

r = shared(np.array([28, 8, -3, 7, -1, 1, 18, 12, 15], dtype=np.float32))
sigma1 = shared(np.array(.2, dtype=np.float32))
alpha0 = shared(np.array(.5, dtype=np.float32))

with Model() as garch:
    alpha1 = Normal('alpha1', 0., 1., dtype='float32')
    BoundedNormal = Bound(Normal, lower=0., upper=(1 - alpha1))
    beta1 = BoundedNormal('beta1', 0., sd=1e3, dtype='float32')

    obs = GARCH11('garchy_garch', omega=alpha0, alpha_1=alpha1,
                  beta_1=beta1, initial_vol=sigma1, observed=r,
                  dtype='float32')
    start = find_MAP(fmin=optimize.fmin_bfgs)

def run(n=1000):
    if n == "short":
        n = 50
    with garch:
        h = find_hessian(start)
        step = Metropolis(garch.vars, h, blocked=True)
        trace = sample(n, step, start)


    burn = n / 10

    traceplot(trace[burn:])
    plots.summary(trace[burn:])


if __name__ == '__main__':
    run()

How are you getting a traceplot or whatever?

[springcoil]

Turns out this works.

from pymc3 import Normal, sample, Model, Bound, traceplot, plots, find_MAP
from pymc3 import *
from pymc3.distributions.timeseries import GARCH11
from scipy import optimize
import theano.tensor as T
from theano import shared
import numpy as np

"""

Example from STAN - slightly altered

GARCH(1,1) example
It is interesting to note just how much more compact this is that the original STAN example

The original implementation is in the STAN documentation by Gelman et al and is reproduced below
A good reference studying GARCH models in depth is  http://www.stat-d.si/mz/mz2.1/posedel.pdf
data {
  int<lower=0> T;
  real r[T];
  real<lower=0> sigma1;
}
parameters {
  real mu;
  real<lower=0> alpha0;
  real<lower=0,upper=1> alpha1;
  real<lower=0, upper=(1-alpha1)> beta1;
}
transformed parameters {
  real<lower=0> sigma[T];
  sigma[1] <- sigma1;
  for (t in 2:T)
    sigma[t] <- sqrt(alpha0
                     + alpha1 * pow(r[t-1] - mu, 2)
                     + beta1 * pow(sigma[t-1], 2));
}
model {
  r ~ normal(mu,sigma);
}
Ported to PyMC3 by Peadar Coyle and Olaf Smits (c) 2016.
"""
r = shared(np.array([28, 8, -3, 7, -1, 1, 18, 12, 15], dtype=np.float32))
sigma1 = shared(np.array(.2, dtype=np.float32))
alpha0 = shared(np.array(.5, dtype=np.float32))

with Model() as garch:
    alpha1 = Normal('alpha1', 0., 1., dtype='float32')
    BoundedNormal = Bound(Normal, lower=0., upper=(1 - alpha1))
    beta1 = BoundedNormal('beta1', 0., sd=1e3, dtype='float32')

    obs = GARCH11('garchy_garch', omega=alpha0, alpha_1=alpha1,
                  beta_1=beta1, initial_vol=sigma1, observed=r,
                  dtype='float32')



def run(n=1000):
    with garch:
        start = find_MAP(fmin=optimize.fmin_bfgs)
        trace = sample(n, Slice(), start=start)


    traceplot(trace)
    plots.summary(trace)


if __name__ == '__main__':
    run()

[olafSmits]

The thing that is missing here is an easy way to obtain the volatility time series for each sample. For each sample the vol series is constructed, but then discarded. Is it possible to save these result somehow? If not, then maybe a function that can generate the vol series for each sample?

[springcoil]

I'm not sure how to answer your question olaf. I'm sure there is a way but I don't know now. If someone can help me debug the model error I added as a line note - then I think this could go into the code base. It'd be the 'minimum' test that this works.

[springcoil]

Ok I'll merge in this commit - and then fix up my other pull request.