fix GARCH11: missing [] in tt.concatenate()

pymc3/distributions/timeseries.py

@@ -186,21 +186,22 @@ class GARCH11(distribution.Continuous):
 
     .. math::
         y_t = \sigma_t * z_t
+
+    .. math::
         \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
+    omega : tensor
+        omega > 0, mean variance
+    alpha_1 : tensor
+        alpha_1 >= 0, autoregressive term coefficient
+    beta_1 : tensor
+        beta_1 >= 0, alpha_1 + beta_1 < 1, moving average term coefficient
+    initial_vol : tensor
+        initial_vol >= 0, initial volatility, sigma_0
     """
 
     def __init__(self, omega, alpha_1, beta_1,
@@ -210,7 +211,7 @@ def __init__(self, omega, alpha_1, beta_1,
         self.omega = omega = tt.as_tensor_variable(omega)
         self.alpha_1 = alpha_1 = tt.as_tensor_variable(alpha_1)
         self.beta_1 = beta_1 = tt.as_tensor_variable(beta_1)
-        self.initial_vol = initial_vol
+        self.initial_vol = tt.as_tensor_variable(initial_vol)
         self.mean = tt.as_tensor_variable(0.)
 
     def get_volatility(self, x):
@@ -224,7 +225,7 @@ def volatility_update(x, vol, w, a, b):
                       outputs_info=[self.initial_vol],
                       non_sequences=[self.omega, self.alpha_1,
                                      self.beta_1])
-        return tt.concatenate(self.initial_vol, vol)
+        return tt.concatenate([[self.initial_vol], vol])
 
     def logp(self, x):
         vol = self.get_volatility(x)

pymc3/tests/test_distributions_timeseries.py

@@ -2,7 +2,7 @@
 
 from ..model import Model
 from ..distributions.continuous import Flat, Normal
-from ..distributions.timeseries import EulerMaruyama, AR
+from ..distributions.timeseries import EulerMaruyama, AR, GARCH11
 from ..sampling import sample, sample_ppc
 from ..theanof import floatX
 
@@ -38,6 +38,31 @@ def test_AR():
 
 
 
+def test_GARCH11():
+    # test data ~ N(0, 1)
+    data = np.array([-1.35078362, -0.81254164,  0.28918551, -2.87043544, -0.94353337,
+                     0.83660719, -0.23336562, -0.58586298, -1.36856736, -1.60832975,
+                     -1.31403141,  0.05446936, -0.97213128, -0.18928725,  1.62011258,
+                     -0.95978616, -2.06536047,  0.6556103 , -0.27816645, -1.26413397])
+    omega = 0.6
+    alpha_1 = 0.4
+    beta_1 = 0.5
+    initial_vol = np.float64(0.9)
+    vol = np.empty_like(data)
+    vol[0] = initial_vol
+    for i in range(len(data)-1):
+        vol[i+1] = np.sqrt(omega + beta_1*vol[i]**2 + alpha_1*data[i]**2)
+
+    with Model() as t:
+        y = GARCH11('y', omega=omega, alpha_1=alpha_1, beta_1=beta_1,
+                    initial_vol=initial_vol, shape=data.shape)
+        z = Normal('z', mu=0, sd=vol, shape=data.shape)
+    garch_like = t['y'].logp({'z':data, 'y': data})
+    reg_like = t['z'].logp({'z':data, 'y': data})
+    np.testing.assert_allclose(garch_like, reg_like)
+
+
+
 def _gen_sde_path(sde, pars, dt, n, x0):
     xs = [x0]
     wt = np.random.normal(size=(n,) if isinstance(x0, float) else (n, x0.size))