Multivariate auto-regressive modeling.
import numpy as np
import matplotlib.pyplot as pp
import nitime.algorithms as alg
import nitime.utils as utils
np.random.seed(1981)
# simulate two mAR systems:
cov = np.diag([0.3, 1.0, 0.2])
# X(t) = 0.8X(t-1) - 0.5X(t-2) + 0.4Z(t-1) + Ex(t)
# Y(t) = 0.9Y(t-1) - 0.8Y(t-2) + Ey(t)
# Z(t) = 0.5Z(t-1) - 0.2Z(t-2) + 0.5Y(t-1) + Ez(t)
a1 = -np.array([[0.8, 0.0, 0.4],
[0.0, 0.9, 0.0],
[0.0, 0.5, 0.5]])
a2 = -np.array([[-0.5, 0.0, 0.0],
[0.0, -0.8, 0.0],
[0.0, 0.0, -0.2]])
a = np.array([a1.copy(), a2.copy()])
# X(t) = 0.8X(t-1) - 0.5X(t-2) + 0.4Z(t-1) + 0.2Y(t-2) + Ex(t)
# Y(t) = 0.9Y(t-1) - 0.8Y(t-2) + Ey(t)
# Z(t) = 0.5Z(t-1) -0.2Z(t-2) + 0.5Y(t-1) + Ez(t)
# just add some feedback from Y to X at 2 lags
a2[0, 1] = -0.2
b = np.array([a1.copy(), a2.copy()])
def extract_ij(i, j, m):
m_ij_rows = m[[i, j]]
return m_ij_rows[:, [i, j]]
w, Haw = alg.transfer_function_xy(a)
w, Hbw = alg.transfer_function_xy(b)
# generate 500 sets of 100 points
N = 500
L = 100
za = np.empty((N, 3, L))
zb = np.empty((N, 3, L))
ea = np.empty((N, 3, L))
eb = np.empty((N, 3, L))
for i in xrange(N):
za[i], ea[i] = utils.generate_mar(a, cov, L)
zb[i], eb[i] = utils.generate_mar(b, cov, L)
# try to estimate the 2nd order (m)AR coefficients--
# average together N estimates of auto-covariance at lags k=0,1,2
Raxx = np.empty((N, 3, 3, 3))
Rbxx = np.empty((N, 3, 3, 3))
for i in xrange(N):
Raxx[i] = utils.autocov_vector(za[i], nlags=3)
Rbxx[i] = utils.autocov_vector(zb[i], nlags=3)
# average trials together to find autocovariance estimate,
# and extract pairwise components from the ac sequence
Raxx = Raxx.mean(axis=0)
xyRa = extract_ij(0, 1, Raxx)
xzRa = extract_ij(0, 2, Raxx)
yzRa = extract_ij(1, 2, Raxx)
Rbxx = Rbxx.mean(axis=0)
xyRb = extract_ij(0, 1, Rbxx)
xzRb = extract_ij(0, 2, Rbxx)
yzRb = extract_ij(1, 2, Rbxx)
# now estimate mAR coefficients and covariance from the full and
# pairwise relationships
Raxx = Raxx.transpose(2, 0, 1)
a_est, cov_est1 = alg.lwr_recursion(Raxx)
a_xy_est, cov_xy_est1 = alg.lwr_recursion(xyRa.transpose(2, 0, 1))
a_xz_est, cov_xz_est1 = alg.lwr_recursion(xzRa.transpose(2, 0, 1))
a_yz_est, cov_yz_est1 = alg.lwr_recursion(yzRa.transpose(2, 0, 1))
Rbxx = Rbxx.transpose(2, 0, 1)
b_est, cov_est2 = alg.lwr_recursion(Rbxx)
b_xy_est, cov_xy_est2 = alg.lwr_recursion(xyRb.transpose(2, 0, 1))
b_xz_est, cov_xz_est2 = alg.lwr_recursion(xzRb.transpose(2, 0, 1))
b_yz_est, cov_yz_est2 = alg.lwr_recursion(yzRb.transpose(2, 0, 1))
fig = pp.figure()
## w, x2y_a, y2x_a = pairwise_causality(0, 1, a_est, cov_est1)
## w, x2y_b, y2x_b = pairwise_causality(0, 1, b_est, cov_est2)
w, x2y_a, y2x_a, _, _ = alg.granger_causality_xy(a_xy_est, cov_xy_est1)
w, x2y_b, y2x_b, _, _ = alg.granger_causality_xy(b_xy_est, cov_xy_est2)
ax = fig.add_subplot(321)
ax.plot(w, x2y_a, 'b--')
ax.plot(w, x2y_b, 'b')
ax.set_title('x to y')
ax.set_ylim((0, 6))
ax = fig.add_subplot(322)
ax.plot(w, y2x_a, 'b--')
ax.plot(w, y2x_b, 'b')
ax.set_title('y to x')
ax.set_ylim((0, 6))
## w, y2z_a, z2y_a = pairwise_causality(1, 2, a_est, cov_est1)
## w, y2z_b, z2y_b = pairwise_causality(1, 2, b_est, cov_est2)
w, y2z_a, z2y_a, _, _ = alg.granger_causality_xy(a_yz_est, cov_yz_est1)
w, y2z_b, z2y_b, _, _ = alg.granger_causality_xy(b_yz_est, cov_yz_est2)
ax = fig.add_subplot(323)
ax.plot(w, y2z_a, 'b--')
ax.plot(w, y2z_b, 'b')
ax.set_title('y to z')
ax.set_ylim((0, 6))
ax = fig.add_subplot(324)
ax.plot(w, z2y_a, 'b--')
ax.plot(w, z2y_b, 'b')
ax.set_title('z to y')
ax.set_ylim((0, 6))
## w, x2z_a, z2x_a = pairwise_causality(0, 2, a_est, cov_est1)
## w, x2z_b, z2x_b = pairwise_causality(0, 2, b_est, cov_est2)
w, x2z_a, z2x_a, _, _ = alg.granger_causality_xy(a_xz_est, cov_xz_est1)
w, x2z_b, z2x_b, _, _ = alg.granger_causality_xy(b_xz_est, cov_xz_est2)
ax = fig.add_subplot(325)
ax.plot(w, x2z_a, 'b--')
ax.plot(w, x2z_b, 'b')
ax.set_title('x to z')
ax.set_ylim((0, 6))
ax = fig.add_subplot(326)
ax.plot(w, z2x_a, 'b--')
ax.plot(w, z2x_b, 'b')
ax.set_title('z to x')
ax.set_ylim((0, 6))
pp.show()
Example source code
You can download the full source code of this example. This same script is also included in the Nitime source distribution under the doc/examples/ directory.