Merge pull request #220 from bashtage/improve-autolag

ENH: Speed up autolag in ADF

Author: bashtage

arch/compat/python.py

@@ -9,49 +9,27 @@
     # always writeable
     from StringIO import StringIO
 
-    BytesIO = StringIO
-    import cPickle
-    pickle = cPickle
-    import urllib2
-    import urlparse
 except ImportError:
     import builtins
-    from io import StringIO, BytesIO
+    from io import StringIO
 
     cStringIO = StringIO
-    import pickle as cPickle
-    pickle = cPickle
-    import urllib.request
-    import urllib.parse
-    from urllib.request import HTTPError, urlretrieve
+
 
 PY3 = sys.version_info[0] == 3
 
 if PY3:
     range = range
     long = int
-    string_types = str,
 
     def lmap(*args, **kwargs):
         return list(map(*args, **kwargs))
 else:
-    string_types = basestring,
     range = xrange
     long = long
     lmap = builtins.map
 
 
-def with_metaclass(meta, *bases):
-    """Create a base class with a metaclass."""
-    # This requires a bit of explanation: the basic idea is to make a dummy
-    # metaclass for one level of class instantiation that replaces itself with
-    # the actual metaclass.
-    class metaclass(meta):
-        def __new__(cls, name, this_bases, d):
-            return meta(name, bases, d)
-    return type.__new__(metaclass, 'temporary_class', (), {})
-
-
 def add_metaclass(metaclass):
     """Class decorator for creating a class with a metaclass."""
     def wrapper(cls):
@@ -80,13 +58,6 @@ def iteritems(obj, **kwargs):
     return func(**kwargs)
 
 
-def iterkeys(obj, **kwargs):
-    func = getattr(obj, "iterkeys", None)
-    if not func:
-        func = obj.keys
-    return func(**kwargs)
-
-
 def itervalues(obj, **kwargs):
     func = getattr(obj, "itervalues", None)
     if not func:

arch/unitroot/critical_values/simulation/adf_z_critical_values_simulation_joblib.py

@@ -11,8 +11,7 @@
 
 from statsmodels.tools.parallel import parallel_func
 import datetime
-from numpy import array, savez, percentile, nan
-from numpy import ones, vstack, arange, cumsum, sum, dot, zeros
+from numpy import array, savez, percentile, nan, ones, vstack, arange, cumsum, sum, dot, zeros
 from numpy.random import RandomState
 from numpy.linalg import pinv
 

arch/unitroot/critical_values/simulation/dfgls_critical_values_simulation.py

@@ -8,7 +8,6 @@
 import datetime
 
 import numpy as np
-from numpy import ones, vstack, arange, diff, cumsum, sqrt, sum
 from numpy.linalg import pinv
 from numpy.random import RandomState
 from statsmodels.compat import range
@@ -57,36 +56,36 @@ def dfgsl_simulation(n, trend, b, rng=None):
     nobs = n
     if trend == 'c':
         c = -7.0
-        z = ones((nobs, 1))
+        z = np.ones((nobs, 1))
     else:
         c = -13.5
-        z = vstack((ones(nobs), arange(1, nobs + 1))).T
+        z = np.vstack((np.ones(nobs), np.arange(1, nobs + 1))).T
 
     ct = c / nobs
 
     delta_z = np.copy(z)
     delta_z[1:, :] = delta_z[1:, :] - (1 + ct) * delta_z[:-1, :]
     delta_z_inv = pinv(delta_z)
     y = standard_normal((n + 50, b))
-    y = cumsum(y, axis=0)
+    y = np.cumsum(y, axis=0)
     y = y[50:, :]
     delta_y = y.copy()
     delta_y[1:, :] = delta_y[1:, :] - (1 + ct) * delta_y[:-1, :]
     detrend_coef = delta_z_inv.dot(delta_y)
     y_detrended = y - z.dot(detrend_coef)
 
-    delta_y_detrended = diff(y_detrended, axis=0)
+    delta_y_detrended = np.diff(y_detrended, axis=0)
     rhs = y_detrended[:-1, :]
     lhs = delta_y_detrended
 
-    xpy = sum(rhs * lhs, 0)
-    xpx = sum(rhs ** 2.0, 0)
+    xpy = np.sum(rhs * lhs, 0)
+    xpx = np.sum(rhs ** 2.0, 0)
     gamma = xpy / xpx
     e = lhs - rhs * gamma
-    sigma2 = sum(e ** 2.0, axis=0) / (n - 1)  # DOF correction?
+    sigma2 = np.sum(e ** 2.0, axis=0) / (n - 1)  # DOF correction?
     gamma_var = sigma2 / xpx
 
-    stat = gamma / sqrt(gamma_var)
+    stat = gamma / np.sqrt(gamma_var)
     return stat
 
 

arch/unitroot/unitroot.py

@@ -1,34 +1,32 @@
 from __future__ import absolute_import, division
-from arch.compat.python import add_metaclass, range, lmap, long
 
 import warnings
 
 from numpy import (diff, ceil, power, sqrt, sum, cumsum, int32, int64, interp, pi,
                    array, inf, abs, log, sort, polyval, empty, argwhere, arange, squeeze)
-from numpy.linalg import pinv
-from scipy.stats import norm
+from numpy.linalg import pinv, qr, inv, solve
 from pandas import DataFrame
-
-from statsmodels.regression.linear_model import OLS
-from statsmodels.tsa.tsatools import lagmat
-from statsmodels.tsa.stattools import _autolag
+from scipy.stats import norm
 from statsmodels.iolib.summary import Summary
 from statsmodels.iolib.table import SimpleTable
+from statsmodels.regression.linear_model import OLS
+from statsmodels.tsa.stattools import _autolag
+from statsmodels.tsa.tsatools import lagmat
 
-from arch.utility import cov_nw
-from arch.utility.exceptions import InvalidLengthWarning, invalid_length_doc
+from arch.compat.python import add_metaclass, range, lmap, long
+from arch.unitroot.critical_values.dfgls import (dfgls_large_p, dfgls_small_p,
+                                                 dfgls_tau_max, dfgls_tau_min,
+                                                 dfgls_tau_star, dfgls_cv_approx)
 from arch.unitroot.critical_values.dickey_fuller import (adf_z_cv_approx, adf_z_large_p, adf_z_max,
                                                          adf_z_min, adf_z_small_p, adf_z_star,
                                                          tau_2010, tau_large_p, tau_max,
                                                          tau_min, tau_small_p, tau_star)
 from arch.unitroot.critical_values.kpss import kpss_critical_values
-from arch.unitroot.critical_values.dfgls import (dfgls_large_p, dfgls_small_p,
-                                                 dfgls_tau_max, dfgls_tau_min,
-                                                 dfgls_tau_star, dfgls_cv_approx)
+from arch.utility import cov_nw
 from arch.utility.array import ensure1d, DocStringInheritor
+from arch.utility.exceptions import InvalidLengthWarning, invalid_length_doc
 from arch.utility.timeseries import add_trend
 
-
 __all__ = ['ADF', 'DFGLS', 'PhillipsPerron', 'KPSS', 'VarianceRatio',
            'kpss_crit', 'mackinnoncrit', 'mackinnonp']
 
@@ -78,27 +76,23 @@ def _autolag_ols(endog, exog, startlag, maxlag, method, modargs=(), regresults=F
     if regresults:
         return _autolag(OLS, endog, exog, startlag, maxlag, method, regresults=regresults)
 
-    resid = squeeze(endog.copy())
-    x = exog[:, startlag:].copy()
+    q, r = qr(exog)
+    qpy = q.T.dot(endog)
+    ypy = endog.T.dot(endog)
+    xpx = exog.T.dot(exog)
+
     sigma2 = empty(maxlag + 1)
     tstat = empty(maxlag + 1)
-    if len(exog) > 0 and startlag > 0:
-        _x = exog[:, :startlag]
-        resid -= _x.dot(pinv(_x).dot(resid))
-        x -= _x.dot(pinv(_x).dot(x))
-    sigma2[0] = (resid ** 2).mean()
+    nobs = float(endog.shape[0])
     tstat[0] = inf
+    for i in range(startlag, startlag + maxlag + 1):
+        b = solve(r[:i, :i], qpy[:i])
+        sigma2[i - startlag] = (ypy - b.T.dot(xpx[:i, :i]).dot(b)) / nobs
+        if method == 't-stat' and i > startlag:
+            xpxi = inv(xpx[:i, :i])
+            stderr = sqrt(sigma2[i - startlag] * xpxi[-1, -1])
+            tstat[i - startlag] = b[-1] / stderr
 
-    for i in range(maxlag):
-        _x = x[:, i:i + 1]
-        xpx = _x.T.dot(_x)
-        beta = squeeze(_x.T.dot(resid) / xpx)
-        resid -= squeeze(beta * _x)
-        x[:, i + 1:] -= _x.dot(_x.T.dot(x[:, i + 1:]) / xpx)
-        sigma2[i + 1] = (resid ** 2).mean()
-        tstat[i + 1] = beta / sqrt(sigma2[i + 1] / xpx)
-
-    nobs = float(resid.shape[0])
     llf = -nobs / 2.0 * (log(2 * pi) + log(sigma2) + 1)
 
     if method == 'aic':

arch/univariate/base.py

@@ -9,8 +9,6 @@
 import warnings
 
 import numpy as np
-from numpy.linalg import matrix_rank
-from numpy import ones, zeros, sqrt, diag, empty, ceil
 import scipy.stats as stats
 import pandas as pd
 from statsmodels.tools.decorators import cache_readonly, resettable_cache
@@ -97,7 +95,7 @@ def format_float_fixed(x, max_digits=10, decimal=4):
     if x == 0:
         return ('{:0.' + str(decimal) + 'f}').format(0.0)
     scale = np.log10(np.abs(x))
-    scale = np.sign(scale) * ceil(np.abs(scale))
+    scale = np.sign(scale) * np.ceil(np.abs(scale))
     if scale > (max_digits - 2 - decimal) or scale < -(decimal - 2):
         formatted = (
             '{0:' + str(max_digits) + '.' + str(decimal) + 'e}').format(x)
@@ -123,7 +121,7 @@ def implicit_constant(x):
         the array has a set of columns that adds to a constant value
     """
     nobs = x.shape[0]
-    rank = matrix_rank(np.hstack((ones((nobs, 1)), x)))
+    rank = np.linalg.matrix_rank(np.hstack((np.ones((nobs, 1)), x)))
     return rank == x.shape[1]
 
 
@@ -149,7 +147,7 @@ def __init__(self, y=None, volatility=None, distribution=None,
         if y is not None:
             self._y_series = ensure1d(y, 'y', series=True)
         else:
-            self._y_series = ensure1d(empty((0,)), 'y', series=True)
+            self._y_series = ensure1d(np.empty((0,)), 'y', series=True)
 
         self._y = np.asarray(self._y_series)
         self._y_original = y
@@ -187,7 +185,7 @@ def constraints(self):
         -----
         Parameters satisfy a.dot(parameters) - b >= 0
         """
-        return empty((0, self.num_params)), empty(0)
+        return np.empty((0, self.num_params)), np.empty(0)
 
     def bounds(self):
         """
@@ -444,7 +442,7 @@ def fit(self, update_freq=1, disp='final', starting_values=None,
         sv_volatility = v.starting_values(resids)
         self._var_bounds = var_bounds = v.variance_bounds(resids)
         v.compute_variance(sv_volatility, resids, sigma2, backcast, var_bounds)
-        std_resids = resids / sqrt(sigma2)
+        std_resids = resids / np.sqrt(sigma2)
 
         # 2. Construct constraint matrices from all models and distribution
         constraints = (self.constraints(),
@@ -454,8 +452,8 @@ def fit(self, update_freq=1, disp='final', starting_values=None,
         num_constraints = [c[0].shape[0] for c in constraints]
         num_constraints = np.array(num_constraints)
         num_params = offsets.sum()
-        a = zeros((num_constraints.sum(), num_params))
-        b = zeros(num_constraints.sum())
+        a = np.zeros((num_constraints.sum(), num_params))
+        b = np.zeros(num_constraints.sum())
 
         for i, c in enumerate(constraints):
             r_en = num_constraints[:i + 1].sum()
@@ -1506,7 +1504,7 @@ def std_err(self):
         """
         Parameter standard error
         """
-        return pd.Series(sqrt(diag(self.param_cov)),
+        return pd.Series(np.sqrt(np.diag(self.param_cov)),
                          index=self._names, name='std_err')
 
     @cache_readonly

arch/univariate/mean.py

@@ -8,7 +8,6 @@
 from collections import OrderedDict
 
 import numpy as np
-from numpy import zeros, empty, ones, isscalar, log
 from pandas import DataFrame
 from scipy.optimize import OptimizeResult
 from statsmodels.tools.decorators import cache_readonly
@@ -191,8 +190,8 @@ def _static_gaussian_loglikelihood(resids):
         nobs = resids.shape[0]
         sigma2 = resids.dot(resids) / nobs
 
-        loglikelihood = -0.5 * nobs * log(2 * np.pi)
-        loglikelihood -= 0.5 * nobs * log(sigma2)
+        loglikelihood = -0.5 * nobs * np.log(2 * np.pi)
+        loglikelihood -= 0.5 * nobs * np.log(sigma2)
         loglikelihood -= 0.5 * nobs
 
         return loglikelihood
@@ -346,10 +345,10 @@ def simulate(self, params, nobs, burn=500, initial_value=None, x=None,
         vol = np.sqrt(sim_data[1])
 
         max_lag = np.max(self._lags)
-        y = zeros(nobs + burn)
+        y = np.zeros(nobs + burn)
         if initial_value is None:
             initial_value = 0.0
-        elif not isscalar(initial_value):
+        elif not np.isscalar(initial_value):
             initial_value = ensure1d(initial_value, 'initial_value')
             if initial_value.shape[0] != max_lag:
                 raise ValueError('initial_value has the wrong shape')
@@ -442,7 +441,7 @@ def _reformat_lags(self):
 
             ind = np.lexsort(np.flipud(lags))
             lags = lags[:, ind]
-            test_mat = zeros((lags.shape[1], np.max(lags)))
+            test_mat = np.zeros((lags.shape[1], np.max(lags)))
             for i in range(lags.shape[1]):
                 test_mat[i, lags[0, i]:lags[1, i]] = 1.0
             rank = np.linalg.matrix_rank(test_mat)
@@ -476,24 +475,24 @@ def _init_model(self):
 
         nobs_orig = self._y.shape[0]
         if self.constant:
-            reg_constant = ones((nobs_orig, 1), dtype=np.float64)
+            reg_constant = np.ones((nobs_orig, 1), dtype=np.float64)
         else:
-            reg_constant = ones((nobs_orig, 0), dtype=np.float64)
+            reg_constant = np.ones((nobs_orig, 0), dtype=np.float64)
 
         if self.lags is not None and nobs_orig > 0:
             maxlag = np.max(self.lags)
             lag_array = lagmat(self._y, maxlag)
-            reg_lags = empty((nobs_orig, self._lags.shape[1]),
-                             dtype=np.float64)
+            reg_lags = np.empty((nobs_orig, self._lags.shape[1]),
+                                dtype=np.float64)
             for i, lags in enumerate(self._lags.T):
                 reg_lags[:, i] = np.mean(lag_array[:, lags[0]:lags[1]], 1)
         else:
-            reg_lags = empty((nobs_orig, 0), dtype=np.float64)
+            reg_lags = np.empty((nobs_orig, 0), dtype=np.float64)
 
         if self._x is not None:
             reg_x = self._x
         else:
-            reg_x = empty((nobs_orig, 0), dtype=np.float64)
+            reg_x = np.empty((nobs_orig, 0), dtype=np.float64)
 
         self.regressors = np.hstack((reg_constant, reg_lags, reg_x))
 
@@ -597,7 +596,7 @@ def _fit_no_arch_normal_errors(self, cov_type='robust'):
             param_cov /= nobs
             cov_type = COV_TYPES['classic_ols']
         elif cov_type in ('robust',):
-            scores = zeros((nobs, self.num_params + 1))
+            scores = np.zeros((nobs, self.num_params + 1))
             scores[:, :self.num_params] = x * e[:, None]
             scores[:, -1] = e ** 2.0 - sigma2
             score_cov = scores.T.dot(scores) / nobs