docs: def simulate_linear_model(A, x0, v, ts_length):
      r"""
      This is a separate function for simulating a vector linear system of
      the form
  
      .. math::
  
          x_{t+1} = A x_t + v_t
  
      given :math:`x_0` = x0
  
      Here :math:`x_t` and :math:`v_t` are both n x 1 and :math:`A` is n x n.
  
      The purpose of separating this functionality out is to target it for
      optimization by Numba.  For the same reason, matrix multiplication is
      broken down into for loops.
  
      Parameters
      ----------
      A : array_like or scalar(float)
          Should be n x n
      x0 : array_like
          Should be n x 1.  Initial condition
      v : np.ndarray
          Should be n x ts_length-1.  Its t-th column is used as the time t
          shock :math:`v_t`
      ts_length : int
          The length of the time series
  
      Returns
      --------
      x : np.ndarray
          Time series with ts_length columns, the t-th column being :math:`x_t`
      """
      A = np.asarray(A)
      n = A.shape[0]
      x = np.empty((n, ts_length))
      x[:, 0] = x0
      for t in range(ts_length-1):
          # x[:, t+1] = A.dot(x[:, t]) + v[:, t]
          for i in range(n):
              x[i, t+1] = v[i, t]                   # Shock
              for j in range(n):
                  x[i, t+1] += A[i, j] * x[j, t]   # Dot Product
      return x