algorithms.autoregressive

Module: algorithms.autoregressive

Autoregressive (AR) processes are processes of the form:

x(n) = a(1)x(n-1) + a(2)x(n-2) + ... + a(P)x(n-P) + e(n)

where e(n) is a white noise process. The usage of ‘e’ suggests interpreting the linear combination of P past values of x(n) as the minimum mean square error linear predictor of x(n) Thus

e(n) = x(n) - a(1)x(n-1) - a(2)x(n-2) - ... - a(P)x(n-P)

Due to whiteness, e(n) is also pointwise uncorrelated–ie,

\begin{align*}
   \text{(i)}   && E\{e(n)e^{*}(n-m)\}& = \delta(n-m) &\\
   \text{(ii)}  && E\{e(n)x^{*}(m)\}  & = 0           & m\neq n\\
   \text{(iii)} && E\{|e|^{2}\} = E\{e(n)e^{*}(n)\} &= E\{e(n)x^{*}(n)\} &
\end{align*}

These principles form the basis of the methods in this module for estimating the AR coefficients and the error/innovations power.

Functions

nitime.algorithms.autoregressive.AR_est_LD(x, order, rxx=None)

Levinson-Durbin algorithm for solving the Hermitian Toeplitz system of Yule-Walker equations in the AR estimation problem

T^{(p)}a^{(p)} = \gamma^{(p+1)}

where

\begin{align*}
T^{(p)} &= \begin{pmatrix}
   R_{0} & R_{1}^{*} & \cdots & R_{p-1}^{*}\\
   R_{1} & R_{0} & \cdots & R_{p-2}^{*}\\
   \vdots & \vdots & \ddots & \vdots\\
   R_{p-1}^{*} & R_{p-2}^{*} & \cdots & R_{0}
\end{pmatrix}\\
a^{(p)} &=\begin{pmatrix} a_1 & a_2 & \cdots a_p \end{pmatrix}^{T}\\
\gamma^{(p+1)}&=\begin{pmatrix}R_1 & R_2 & \cdots & R_p \end{pmatrix}^{T}
\end{align*}

and R_k is the autocorrelation of the kth lag

Parameters:

x : ndarray

the zero-mean stochastic process

order : int

the AR model order–IE the rank of the system.

rxx : ndarray, optional

(at least) order+1 samples of the autocorrelation sequence

Returns:

ak, sig_sq

The AR coefficients for 1 <= k <= p, and the variance of the driving white noise process

nitime.algorithms.autoregressive.AR_est_YW(x, order, rxx=None)

Determine the autoregressive (AR) model of a random process x using the Yule Walker equations. The AR model takes this convention:

x(n) = a(1)x(n-1) + a(2)x(n-2) + \dots + a(p)x(n-p) + e(n)

where e(n) is a zero-mean white noise process with variance sig_sq, and p is the order of the AR model. This method returns the a_i and sigma

The orthogonality property of minimum mean square error estimates states that

E\{e(n)x^{*}(n-k)\} = 0 \quad 1\leq k\leq p

Inserting the definition of the error signal into the equations above yields the Yule Walker system of equations:

R_{xx}(k) = \sum_{i=1}^{p}a(i)R_{xx}(k-i) \quad1\leq k\leq p

Similarly, the variance of the error process is

E\{e(n)e^{*}(n)\}   = E\{e(n)x^{*}(n)\} = R_{xx}(0)-\sum_{i=1}^{p}a(i)R^{*}(i)

Parameters:

x : ndarray

The sampled autoregressive random process

order : int

The order p of the AR system

rxx : ndarray (optional)

An optional, possibly unbiased estimate of the autocorrelation of x

Returns:

ak, sig_sq : The estimated AR coefficients and innovations variance

nitime.algorithms.autoregressive.AR_psd(ak, sigma_v, n_freqs=1024, sides='onesided')

Compute the PSD of an AR process, based on the process coefficients and covariance

n_freqs : int
The number of spacings on the frequency grid from [-PI,PI). If sides==’onesided’, n_freqs/2+1 frequencies are computed from [0,PI]
sides : str (optional)
Indicates whether to return a one-sided or two-sided PSD
Returns:

(w, ar_psd)

w : Array of normalized frequences from [-.5, .5) or [0,.5]

ar_psd : A PSD estimate computed by sigma_v / |1-a(f)|**2 , where

a(f) = DTFT(ak)

nitime.algorithms.autoregressive.MAR_est_LWR(x, order, rxx=None)

MAR estimation, using the LWR algorithm, as in Morf et al.

Parameters:

x : ndarray

The sampled autoregressive random process

order : int

The order P of the AR system

rxx : ndarray (optional)

An optional, possibly unbiased estimate of the autocovariance of x

Returns:

a, ecov : The system coefficients and the estimated covariance

nitime.algorithms.autoregressive.coherence_from_spectral(Sw)

Compute the spectral coherence between processes X and Y, given their spectral matrix S(w)

Parameters:

Sw : ndarray

spectral matrix

nitime.algorithms.autoregressive.granger_causality_xy(a, cov, n_freqs=1024)

Compute the Granger causality between processes X and Y, which are linked in a multivariate autoregressive (mAR) model parameterized by coefficient matrices a(i) and the innovations covariance matrix

X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

Parameters:

a : ndarray, (P,2,2)

coefficient matrices characterizing the autoregressive mixing

cov : ndarray, (2,2)

covariance matrix characterizing the innovations vector

n_freqs: int

number of frequencies to compute in the fourier transform

Returns:

w, f_x_on_y, f_y_on_x, f_xy, Sw

  1. vector of frequencies
  2. function of the Granger causality of X on Y
  3. function of the Granger causality of Y on X
  4. function of the ‘instantaneous causality’ between X and Y
  5. spectral density matrix
nitime.algorithms.autoregressive.interdependence_xy(Sw)

Compute the ‘total interdependence’ between processes X and Y, given their spectral matrix S(w)

Parameters:

Sw : ndarray

spectral matrix

Returns:

fxy(w)

interdependence function of frequency

nitime.algorithms.autoregressive.lwr_recursion(r)

Perform a Levinson-Wiggins[Whittle]-Robinson recursion to find the coefficients a(i) that satisfy the matrix version of the Yule-Walker system of P + 1 equations:

sum_{i=0}^{P} a(i)r(k-i) = 0, for k = {1,2,...,P}

with the additional equation

sum_{i=0}^{P} a(i)r(-k) = V

where V is the covariance matrix of the innovations process, and a(0) is fixed at the identity matrix

Also note that r is defined as:

r(k) = E{ X(t)X*(t-k) } ( * = conjugate transpose ) r(-k) = r*(k)

This routine adapts the algorithm found in eqs (1)-(11) in Morf, Vieira, Kailath 1978

Parameters:

r : ndarray, shape (P + 1, nc, nc)

Returns:

a : ndarray (P,nc,nc)

coefficient sequence of order P

sigma : ndarray (nc,nc)

covariance estimate

nitime.algorithms.autoregressive.spectral_matrix_xy(Hw, cov)

Compute the spectral matrix S(w), from the convention:

X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

The formulation follows from Ding, Chen, Bressler 2008, pg 6 eqs (11) to (15)

The transfer function H(w) should be computed first from transfer_function_xy()

Parameters:

Hw : ndarray (2, 2, n_freqs)

Pre-computed transfer function from transfer_function_xy()

cov : ndarray (2, 2)

The covariance between innovations processes in Err[t]

Returns:

Sw : ndarrays

matrix of spectral density functions

nitime.algorithms.autoregressive.transfer_function_xy(a, n_freqs=1024)

Helper routine to compute the transfer function H(w) based on sequence of coefficient matrices A(i). The z transforms follow from this definition:

X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

Parameters:

a : ndarray, shape (P, 2, 2)

sequence of coef matrices describing an mAR process

n_freqs : int, optional

number of frequencies to compute in range [0,PI]

Returns:

Hw : ndarray

The transfer function from innovations process vector to mAR process X