MAINT: Distinguish between private and public namespaces (#601)

* MAINT: Rename `ecdf.py` to `_ecdf.py`

* MAINT: Add `ecdf.py` with warnings to avoid breaking changes

* MAINT: Rename files

* MAINT: Write deprecation warnings

* MAINT: fix errors and imports

* MAINT: Rename files

* MAINT: Use imports and add deprecation warnings

* MAINT: Rename files

* MAINT: Fix imports and add deprecation warnings

* MAINT: Rename files

* MAINT: Fix imports

* Fix merge issues

* MAINT: Correct version to 0.8.0

* Fix docs warning

quantecon/__init__.py

@@ -12,42 +12,51 @@
 		"Cannot import numba from current anaconda distribution. \
 			Please run `conda install numba` to install the latest version.")
 
-#-Modules-#
+# Modules
 from . import distributions
 from . import game_theory
 from . import quad
 from . import random
 from . import optimize
 
-#-Objects-#
-from .compute_fp import compute_fixed_point
-from .discrete_rv import DiscreteRV
-from .dle import DLE
-from .ecdf import ECDF
-from .estspec import smooth, periodogram, ar_periodogram
-# from .game_theory import <objects-here> 							#Place Holder if we wish to promote any general objects to the qe namespace.
-from .graph_tools import DiGraph, random_tournament_graph
-from .gridtools import (
-	cartesian, mlinspace, cartesian_nearest_index, simplex_grid, simplex_index
-)
-from .inequality import lorenz_curve, gini_coefficient, shorrocks_index, \
-	rank_size
-from .kalman import Kalman
-from .lae import LAE
-from .arma import ARMA
-from .lqcontrol import LQ, LQMarkov
-from .filter import hamilton_filter
-from .lqnash import nnash
-from .lss import LinearStateSpace
-from .matrix_eqn import solve_discrete_lyapunov, solve_discrete_riccati
-from .quadsums import var_quadratic_sum, m_quadratic_sum
-#->Propose Delete From Top Level
-#Promote to keep current examples working
+# Objects
+from ._compute_fp import compute_fixed_point
+from ._discrete_rv import DiscreteRV
+from ._dle import DLE
+from ._ecdf import ECDF
+from ._estspec import smooth, periodogram, ar_periodogram
+from ._graph_tools import DiGraph, random_tournament_graph
+from ._gridtools import (cartesian, mlinspace, simplex_grid, simplex_index,
+						 num_compositions)
+from ._inequality import (lorenz_curve, gini_coefficient, shorrocks_index,
+						  rank_size)
+from ._kalman import Kalman
+from ._lae import LAE
+from ._arma import ARMA
+from ._lqcontrol import LQ, LQMarkov
+from ._filter import hamilton_filter
+from ._lqnash import nnash
+from ._ivp import IVP
+from ._lss import LinearStateSpace
+from ._matrix_eqn import solve_discrete_lyapunov, solve_discrete_riccati
+from ._quadsums import var_quadratic_sum, m_quadratic_sum
+from ._rank_nullspace import rank_est, nullspace
+from ._robustlq import RBLQ
+
+#-> Propose Delete From Top Level
+# Promote to keep current examples working
 from .markov import MarkovChain, random_markov_chain, random_stochastic_matrix, \
 	gth_solve, tauchen, rouwenhorst
-#Imports that Should be Deprecated with markov package
-from .markov import mc_compute_stationary, mc_sample_path
-#<-
-from .rank_nullspace import rank_est, nullspace
-from .robustlq import RBLQ
+
+
 from .util import searchsorted, fetch_nb_dependencies, tic, tac, toc
+#<-
+
+# Imports that should be deprecated with markov package
+from .markov import mc_compute_stationary, mc_sample_path
+
+# Imports that are deprecated and will be removed in further versions
+from . import (ecdf, arma, compute_fp, discrete_rv, dle, estspec, filter,
+			   graph_tools, gridtools, inequality, ivp, kalman, lae,
+			   lqcontrol, lqnash, lss, matrix_eqn, quadsums, rank_nullspace,
+			   robustlq)

quantecon/_arma.py

@@ -0,0 +1,258 @@
+"""
+Provides functions for working with and visualizing scalar ARMA processes.
+
+TODO: 1. Fix warnings concerning casting complex variables back to floats
+
+"""
+import numpy as np
+from .util import check_random_state
+
+
+class ARMA:
+    r"""
+    This class represents scalar ARMA(p, q) processes.
+
+    If phi and theta are scalars, then the model is
+    understood to be
+
+    .. math::
+
+        X_t = \phi X_{t-1} + \epsilon_t + \theta \epsilon_{t-1}
+
+    where :math:`\epsilon_t` is a white noise process with standard
+    deviation :math:`\sigma`.  If phi and theta are arrays or sequences,
+    then the interpretation is the ARMA(p, q) model
+
+    .. math::
+
+        X_t = \phi_1 X_{t-1} + ... + \phi_p X_{t-p} +
+
+        \epsilon_t + \theta_1 \epsilon_{t-1} + ...  +
+        \theta_q \epsilon_{t-q}
+
+    where
+
+        * :math:`\phi = (\phi_1, \phi_2,..., \phi_p)`
+        * :math:`\theta = (\theta_1, \theta_2,..., \theta_q)`
+        * :math:`\sigma` is a scalar, the standard deviation of the
+          white noise
+
+    Parameters
+    ----------
+    phi : scalar or iterable or array_like(float)
+        Autocorrelation values for the autocorrelated variable.
+        See above for explanation.
+    theta : scalar or iterable or array_like(float)
+        Autocorrelation values for the white noise of the model.
+        See above for explanation
+    sigma : scalar(float)
+        The standard deviation of the white noise
+
+    Attributes
+    ----------
+    phi, theta, sigma : see Parmeters
+    ar_poly : array_like(float)
+        The polynomial form that is needed by scipy.signal to do the
+        processing we desire.  Corresponds with the phi values
+    ma_poly : array_like(float)
+        The polynomial form that is needed by scipy.signal to do the
+        processing we desire.  Corresponds with the theta values
+
+    """
+    def __init__(self, phi, theta=0, sigma=1):
+        self._phi, self._theta = phi, theta
+        self.sigma = sigma
+        self.set_params()
+
+    def __repr__(self):
+        m = "ARMA(phi=%s, theta=%s, sigma=%s)"
+        return m % (self.phi, self.theta, self.sigma)
+
+    def __str__(self):
+        m = "An ARMA({p}, {q}) process"
+        p = np.asarray(self.phi).size
+        q = np.asarray(self.theta).size
+        return m.format(p=p, q=q)
+
+    # Special latex print method for working in notebook
+    def _repr_latex_(self):
+        m = r"$X_t = "
+        phi = np.atleast_1d(self.phi)
+        theta = np.atleast_1d(self.theta)
+        rhs = ""
+        for (tm, phi_p) in enumerate(phi):
+            # don't include terms if they are equal to zero
+            if abs(phi_p) > 1e-12:
+                rhs += r"%+g X_{t-%i}" % (phi_p, tm+1)
+
+        if rhs[0] == "+":
+            rhs = rhs[1:]  # remove initial `+` if phi_1 was positive
+
+        rhs += r" + \epsilon_t"
+
+        for (tm, th_q) in enumerate(theta):
+            # don't include terms if they are equal to zero
+            if abs(th_q) > 1e-12:
+                rhs += r"%+g \epsilon_{t-%i}" % (th_q, tm+1)
+
+        return m + rhs + "$"
+
+    @property
+    def phi(self):
+        return self._phi
+
+    @phi.setter
+    def phi(self, new_value):
+        self._phi = new_value
+        self.set_params()
+
+    @property
+    def theta(self):
+        return self._theta
+
+    @theta.setter
+    def theta(self, new_value):
+        self._theta = new_value
+        self.set_params()
+
+    def set_params(self):
+        r"""
+        Internally, scipy.signal works with systems of the form
+
+        .. math::
+
+            ar_{poly}(L) X_t = ma_{poly}(L) \epsilon_t
+
+        where L is the lag operator. To match this, we set
+
+        .. math::
+
+            ar_{poly} = (1, -\phi_1, -\phi_2,..., -\phi_p)
+
+            ma_{poly} = (1, \theta_1, \theta_2,..., \theta_q)
+
+        In addition, ar_poly must be at least as long as ma_poly.
+        This can be achieved by padding it out with zeros when required.
+
+        """
+        # === set up ma_poly === #
+        ma_poly = np.asarray(self._theta)
+        self.ma_poly = np.insert(ma_poly, 0, 1)  # The array (1, theta)
+
+        # === set up ar_poly === #
+        if np.isscalar(self._phi):
+            ar_poly = np.array(-self._phi)
+        else:
+            ar_poly = -np.asarray(self._phi)
+        self.ar_poly = np.insert(ar_poly, 0, 1)  # The array (1, -phi)
+
+        # === pad ar_poly with zeros if required === #
+        if len(self.ar_poly) < len(self.ma_poly):
+            temp = np.zeros(len(self.ma_poly) - len(self.ar_poly))
+            self.ar_poly = np.hstack((self.ar_poly, temp))
+
+    def impulse_response(self, impulse_length=30):
+        """
+        Get the impulse response corresponding to our model.
+
+        Returns
+        -------
+        psi : array_like(float)
+            psi[j] is the response at lag j of the impulse response.
+            We take psi[0] as unity.
+
+        """
+        from scipy.signal import dimpulse
+        sys = self.ma_poly, self.ar_poly, 1
+        times, psi = dimpulse(sys, n=impulse_length)
+        psi = psi[0].flatten()  # Simplify return value into flat array
+
+        return psi
+
+    def spectral_density(self, two_pi=True, res=1200):
+        r"""
+        Compute the spectral density function.  The spectral density is
+        the discrete time Fourier transform of the autocovariance
+        function.  In particular,
+
+        .. math::
+
+            f(w) = \sum_k \gamma(k) \exp(-ikw)
+
+        where gamma is the autocovariance function and the sum is over
+        the set of all integers.
+
+        Parameters
+        ----------
+        two_pi : Boolean, optional
+            Compute the spectral density function over :math:`[0, \pi]` if
+            two_pi is False and :math:`[0, 2 \pi]` otherwise.  Default value is
+            True
+        res : scalar or array_like(int), optional(default=1200)
+            If res is a scalar then the spectral density is computed at
+            `res` frequencies evenly spaced around the unit circle, but
+            if res is an array then the function computes the response
+            at the frequencies given by the array
+
+        Returns
+        -------
+        w : array_like(float)
+            The normalized frequencies at which h was computed, in
+            radians/sample
+        spect : array_like(float)
+            The frequency response
+
+        """
+        from scipy.signal import freqz
+        w, h = freqz(self.ma_poly, self.ar_poly, worN=res, whole=two_pi)
+        spect = h * np.conj(h) * self.sigma**2
+
+        return w, spect
+
+    def autocovariance(self, num_autocov=16):
+        """
+        Compute the autocovariance function from the ARMA parameters
+        over the integers range(num_autocov) using the spectral density
+        and the inverse Fourier transform.
+
+        Parameters
+        ----------
+        num_autocov : scalar(int), optional(default=16)
+            The number of autocovariances to calculate
+
+        """
+        spect = self.spectral_density()[1]
+        acov = np.fft.ifft(spect).real
+
+        # num_autocov should be <= len(acov) / 2
+        return acov[:num_autocov]
+
+    def simulation(self, ts_length=90, random_state=None):
+        """
+        Compute a simulated sample path assuming Gaussian shocks.
+
+        Parameters
+        ----------
+        ts_length : scalar(int), optional(default=90)
+            Number of periods to simulate for
+
+        random_state : int or np.random.RandomState/Generator, optional
+            Random seed (integer) or np.random.RandomState or Generator
+            instance to set the initial state of the random number
+            generator for reproducibility. If None, a randomly
+            initialized RandomState is used.
+
+        Returns
+        -------
+        vals : array_like(float)
+            A simulation of the model that corresponds to this class
+
+        """
+        from scipy.signal import dlsim
+        random_state = check_random_state(random_state)
+
+        sys = self.ma_poly, self.ar_poly, 1
+        u = random_state.standard_normal((ts_length, 1)) * self.sigma
+        vals = dlsim(sys, u)[1]
+
+        return vals.flatten()