Skip to content

Add GARCH11 model #987

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Feb 21, 2016
Merged

Add GARCH11 model #987

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 51 additions & 1 deletion pymc3/distributions/timeseries.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
import theano.tensor as T
from theano import scan

from .continuous import Normal, Flat
from .distribution import Continuous

__all__ = ['AR1', 'GaussianRandomWalk']
__all__ = ['AR1', 'GaussianRandomWalk', 'GARCH11']


class AR1(Continuous):
Expand Down Expand Up @@ -71,3 +72,52 @@ def logp(self, x):

innov_like = Normal.dist(mu=x_im1 + mu, tau=tau, sd=sd).logp(x_i)
return init.logp(x[0]) + T.sum(innov_like)


class GARCH11(Continuous):
"""
GARCH(1,1) with Normal innovations. The model is specified by

y_t = sigma_t * z_t
sigma_t^2 = omega + alpha_1 * y_{t-1}^2 + beta_1 * sigma_{t-1}^2

with z_t iid and Normal with mean zero and unit standard deviation.

Parameters
----------
omega : distribution
omega > 0, distribution for mean variance
alpha_1 : distribution
alpha_1 >= 0, distribution for autoregressive term
beta_1 : distribution
beta_1 >= 0, alpha_1 + beta_1 < 1, distribution for moving
average term
initial_vol : distribution
initial_vol >= 0, distribution for initial volatility, sigma_0
"""
def __init__(self, omega=None, alpha_1=None, beta_1=None,
initial_vol=None, *args, **kwargs):
super(GARCH11, self).__init__(*args, **kwargs)

self.omega = omega
self.alpha_1 = alpha_1
self.beta_1 = beta_1
self.initial_vol = initial_vol
self.mean = 0

def _get_volatility(self, x):

def volatility_update(x, vol, w, a, b):
return T.sqrt(w + a * T.square(x) + b * T.square(vol))

vol, _ = scan(fn=volatility_update,
sequences=[x],
outputs_info=[self.initial_vol],
non_sequences=[self.omega, self.alpha_1,
self.beta_1])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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'>)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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().

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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()

return vol

def logp(self, x):
vol = self._get_volatility(x[:-1])
return (Normal.dist(0., sd=self.initial_vol).logp(x[0]) +
T.sum(Normal.dist(0, sd=vol).logp(x[1:])))