{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time-domain Granger Causality\n",
    "## Pairwise Granger Causality\n",
    "The Granger causality is a method to determine functional connectivity between time-series using autoregressive modelling. In the simpliest pairwise Granger causality case for signals X and Y the data are modelled as autoregressive processes. Each of these processes has two representations. The first representation contains the history of the signal X itself and a prediction error (or noise a.k.a. residual), whereas the second also incorporates the history of the other signal. \n",
    "\n",
    "If inclusion of the history of Y next to the history of X into X model reduces the prediction error compared to just the history of X alone, Y is said to Granger cause X. The same can be done by interchanging the signals to determine if X Granger causes Y.\n",
    "\n",
    "## Conditional Granger Causality\n",
    "Conditional Granger causality can be used to further investigate this functional connectivity. Given signals X, Y and Z, we find that Y Granger causes X, but we want to test if this causality is mediated through Z. We can use Z as a condition for the aforementioned Granger causality.\n",
    "\n",
    "In order to illustrate the function of time-domain Granger causality we will be using examples from Ding et al. (2006) chapter. Specifically, we will have two cases of three signals. In the first case we will have indirect connectivity only, whereas in the second case both direct and indirect connectivities will be present.\n",
    "\n",
    "References: Ding M., Chen Y. and Bressler S.L. (2006) Granger Causality: Basic Theory and Application to Neuroscience. https://arxiv.org/abs/q-bio/0608035\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from elephant.causality.granger import pairwise_granger, conditional_granger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)\n",
    "# Indirect causal influence diagram\n",
    "node1 = plt.Circle((0.2, 0.2), 0.1, color='red')\n",
    "node2 = plt.Circle((0.5, 0.6), 0.1, color='red')\n",
    "node3 = plt.Circle((0.8, 0.2), 0.1, color='red')\n",
    "ax1.set_aspect(1)\n",
    "ax1.arrow(0.28, 0.3, 0.1, 0.125, width=0.02, color='k')\n",
    "ax1.arrow(0.6, 0.5, 0.1, -0.125, width=0.02, color='k')\n",
    "ax1.add_artist(node1)\n",
    "ax1.add_artist(node2)\n",
    "ax1.add_artist(node3)\n",
    "ax1.text(0.2, 0.2, 'Y', horizontalalignment='center', verticalalignment='center')\n",
    "ax1.text(0.5, 0.6, 'Z', horizontalalignment='center', verticalalignment='center')\n",
    "ax1.text(0.8, 0.2, 'X', horizontalalignment='center', verticalalignment='center')\n",
    "ax1.set_title('Indirect only')\n",
    "ax1.set_xbound((0, 1))\n",
    "ax1.set_ybound((0, 0.8))\n",
    "\n",
    "ax1.tick_params(axis='both', which='both', bottom=False, top=False, labelbottom=False, \n",
    "                right=False, left=False, labelleft=False)\n",
    "\n",
    "# Both direct and indirect causal influence diagram\n",
    "node1 = plt.Circle((0.2, 0.2), 0.1, color='g')\n",
    "node2 = plt.Circle((0.5, 0.6), 0.1, color='g')\n",
    "node3 = plt.Circle((0.8, 0.2), 0.1, color='g')\n",
    "ax2.set_aspect(1)\n",
    "ax2.arrow(0.28, 0.3, 0.1, 0.125, width=0.02, color='k')\n",
    "ax2.arrow(0.35, 0.2, 0.2, 0.0, width=0.02, color='k')\n",
    "ax2.arrow(0.6, 0.5, 0.1, -0.125, width=0.02, color='k')\n",
    "ax2.add_artist(node1)\n",
    "ax2.add_artist(node2)\n",
    "ax2.add_artist(node3)\n",
    "ax2.text(0.2, 0.2, 'Y', horizontalalignment='center', verticalalignment='center')\n",
    "ax2.text(0.5, 0.6, 'Z', horizontalalignment='center', verticalalignment='center')\n",
    "ax2.text(0.8, 0.2, 'X', horizontalalignment='center', verticalalignment='center')\n",
    "ax2.set_xbound((0, 1))\n",
    "ax2.set_ybound((0, 0.8))\n",
    "\n",
    "ax2.tick_params(axis='both', which='both', bottom=False, top=False, labelbottom=False, \n",
    "                right=False, left=False, labelleft=False)\n",
    "ax2.set_title('Both direct and indirect')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(length=30000, causality_type=\"indirect\"):\n",
    "    \"\"\"\n",
    "    Recreated from Example 2 section 5.2 of :cite:'granger-Ding06-0608035'.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    length : int\n",
    "        The length of the signals to be generated (i.e. shape(signal) = (3, length))\n",
    "    causality_type: str\n",
    "        Type of causal influence in the data:\n",
    "            'indirect' for indirect causal influence only (i.e. Y -> Z -> X)\n",
    "            'both' for direct and indirect causal influence\n",
    "    \n",
    "    \n",
    "    Notes\n",
    "    -----\n",
    "    Taken from elephant.test.test_causality.ConditionalGrangerTestCase\n",
    "\n",
    "    \"\"\"\n",
    "    if causality_type == \"indirect\":\n",
    "        y_t_lag_2 = 0\n",
    "    elif causality_type == \"both\":\n",
    "        y_t_lag_2 = 0.2\n",
    "    else:\n",
    "        raise ValueError(\"causality_type should be either 'indirect' or \"\n",
    "                         \"'both'\")\n",
    "\n",
    "    order = 2\n",
    "    signal = np.zeros((3, length + order))\n",
    "\n",
    "    weights_1 = np.array([[0.8, 0, 0.4],\n",
    "                          [0, 0.9, 0],\n",
    "                          [0., 0.5, 0.5]])\n",
    "\n",
    "    weights_2 = np.array([[-0.5, y_t_lag_2, 0.],\n",
    "                          [0., -0.8, 0],\n",
    "                          [0, 0, -0.2]])\n",
    "\n",
    "    weights = np.stack((weights_1, weights_2))\n",
    "\n",
    "    noise_covariance = np.array([[0.3, 0.0, 0.0],\n",
    "                                 [0.0, 1., 0.0],\n",
    "                                 [0.0, 0.0, 0.2]])\n",
    "\n",
    "    for i in range(length):\n",
    "        for lag in range(order):\n",
    "            signal[:, i + order] += np.dot(weights[lag],\n",
    "                                           signal[:, i + 1 - lag])\n",
    "        rnd_var = np.random.multivariate_normal([0, 0, 0],\n",
    "                                                noise_covariance)\n",
    "        signal[:, i + order] += rnd_var\n",
    "\n",
    "    signal = signal[:, 2:]\n",
    "\n",
    "    return signal.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Indirect causality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "\n",
    "# Indirect causality\n",
    "xyz_indirect_sig = generate_data(length=10000, causality_type='indirect')\n",
    "xy_indirect_sig = xyz_indirect_sig[:, :2]\n",
    "indirect_pairwise_gc = pairwise_granger(xy_indirect_sig, max_order=10, information_criterion='aic')\n",
    "print(indirect_pairwise_gc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indirect causality (conditioned on z)\n",
    "indirect_cond_gc = conditional_granger(xyz_indirect_sig, max_order=10, information_criterion='aic')\n",
    "print(indirect_cond_gc)\n",
    "print('Zero value indicates total dependence on signal Z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Both direct and indirect causality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Both direct and indirect causality\n",
    "xyz_both_sig = generate_data(length=10000, causality_type='both')\n",
    "xy_both_sig = xyz_both_sig[:, :2]\n",
    "both_pairwise_gc = pairwise_granger(xy_both_sig, max_order=10, information_criterion='aic')\n",
    "print(both_pairwise_gc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Both direct and indirect causality (conditioned on z)\n",
    "both_cond_gc = conditional_granger(xyz_both_sig, max_order=10, information_criterion='aic')\n",
    "print(both_cond_gc)\n",
    "print('Non-zero value indicates the presence of direct Y to X influence')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}