FIX: delayed ARMA impulse response and simulation when p > q

quantecon/_arma.py

@@ -131,8 +131,9 @@ def set_params(self):
 
             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.
+        In addition, ar_poly and ma_poly must have the same length,
+        otherwise scipy.signal delays the output by the difference in
+        length. This is achieved by padding the shorter one with zeros.
 
         """
         # === set up ma_poly === #
@@ -146,10 +147,13 @@ def set_params(self):
             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 === #
+        # === pad the shorter polynomial with zeros === #
         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))
+        elif len(self.ma_poly) < len(self.ar_poly):
+            temp = np.zeros(len(self.ar_poly) - len(self.ma_poly))
+            self.ma_poly = np.hstack((self.ma_poly, temp))
 
     def impulse_response(self, impulse_length=30):
         """

quantecon/tests/test_arma.py

@@ -4,7 +4,7 @@
 
 """
 import numpy as np
-from numpy.testing import assert_array_equal, assert_
+from numpy.testing import assert_allclose, assert_array_equal, assert_
 from quantecon import ARMA
 
 
@@ -42,3 +42,42 @@ def test_impulse_response(self):
         imp_resp = lp.impulse_response(impulse_length=75)
 
         assert_(imp_resp.size == 75)
+
+
+def _ma_coefficients(phi, theta, n):
+    # psi_j = theta_j + sum_i phi_i psi_{j-i}, with psi_0 = 1
+    phi, theta = np.atleast_1d(phi), np.atleast_1d(theta)
+    psi = np.zeros(n)
+    for j in range(n):
+        psi[j] = 1.0 if j == 0 else (theta[j - 1] if j <= len(theta) else 0.0)
+        for i in range(1, min(j, len(phi)) + 1):
+            psi[j] += phi[i - 1] * psi[j - i]
+    return psi
+
+
+ORDERS = [
+    (0.9, 0),                          # AR(1)
+    ([0.5, -0.2], 0),                  # AR(2), default theta
+    ([0.5, 0.1, 0.05], 0),             # AR(3)
+    ([0.5, -0.2], [0.3]),              # ARMA(2, 1)
+    ([0.5], [0.3, 0.2]),               # ARMA(1, 2)
+    ([0.95, -0.4, -0.4], np.zeros(3)),
+]
+
+
+def test_impulse_response_matches_ma_representation():
+    for phi, theta in ORDERS:
+        imp_resp = ARMA(phi, theta).impulse_response(impulse_length=20)
+        assert_allclose(imp_resp, _ma_coefficients(phi, theta, 20))
+
+
+def test_simulation_matches_recursion():
+    for phi, theta in ORDERS:
+        sigma, seed, n = 0.5, 1234, 50
+        sim = ARMA(phi, theta, sigma).simulation(ts_length=n,
+                                                 random_state=seed)
+        u = np.random.RandomState(seed).standard_normal(n) * sigma
+        # With zero initial conditions, X_t = sum_j psi_j u_{t-j}
+        psi = _ma_coefficients(phi, theta, n)
+        expected = np.array([psi[:t + 1] @ u[t::-1] for t in range(n)])
+        assert_allclose(sim, expected)