{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from sklearn import linear_model\n",
    "from sklearn import metrics\n",
    "import scipy\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1000.0, -106.84462699261802, 93.79677127071842)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 575876 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def sim_random_walk(n):\n",
    "    eps = np.random.normal(size=n)\n",
    "    return np.cumsum(eps)\n",
    "\n",
    "n = 1000\n",
    "sns.set(style='white')\n",
    "fig = plt.figure(figsize=(16, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "plt.xlim(0, n)\n",
    "for _ in range(100):\n",
    "    ax.plot(sim_random_walk(n))\n",
    "    \n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7ff87adb1780>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 100324 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(18)\n",
    "n = 500\n",
    "fig, ax = plt.subplots(1,2, figsize=(16, 6), gridspec_kw={'width_ratios':[3, 1]})\n",
    "eps = np.random.normal(size=n)\n",
    "ax[0].plot(eps)\n",
    "sns.distplot(eps, ax=ax[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2],\n",
       "       [1, 2, 3],\n",
       "       [2, 3, 4],\n",
       "       [3, 4, 5],\n",
       "       [4, 5, 6],\n",
       "       [5, 6, 7],\n",
       "       [6, 7, 8]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def lag_view(x, order):\n",
    "    \"\"\"\n",
    "    For every value X_i create a row that lags k values: [X_i-1, X_i-2, ... X_i-k]\n",
    "    \"\"\"\n",
    "    y = x.copy()\n",
    "    # Create features by shifting the window of `order` size by one step.\n",
    "    # This results in a 2D array [[t1, t2, t3], [t2, t3, t4], ... [t_k-2, t_k-1, t_k]]\n",
    "    x = np.array([y[-(i + order):][:order] for i in range(y.shape[0])])\n",
    "    \n",
    "    # Reverse the array as we started at the end and remove duplicates.\n",
    "    # Note that we truncate the features [order -1:] and the labels [order]\n",
    "    # This is the shifting of the features with one time step compared to the labels\n",
    "    x = np.stack(x)[::-1][order - 1: -1]\n",
    "    y = y[order:]\n",
    "\n",
    "    return x, y\n",
    "\n",
    "lag_view(np.arange(10), 3)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 231660 characters]\n",
      "text/plain": [
       "<Figure size 1296x864 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def ma_process(eps, theta):\n",
    "    \"\"\"\n",
    "    Creates an MA(q) process with a zero mean (mean not included in implementation).\n",
    "    :param eps: (array) White noise signal.\n",
    "    :param theta: (array/ list) Parameters of the process.\n",
    "    \"\"\"\n",
    "    # reverse the order of theta as Xt, Xt-1, Xt-k in an array is Xt-k, Xt-1, Xt.\n",
    "    theta = np.array([1] + list(theta))[::-1][:, None]\n",
    "    eps_q, _ = lag_view(eps, len(theta))\n",
    "    return eps_q @ theta\n",
    "\n",
    "fig = plt.figure(figsize=(18, 4 * 3))\n",
    "a = 310\n",
    "for i in range(0, 11, 5):\n",
    "    a += 1\n",
    "    theta = np.random.uniform(0, 1, size=i + 1)\n",
    "    plt.subplot(a)\n",
    "    plt.title(f'$\\\\theta$ = {theta.round(2)}')\n",
    "    plt.plot(ma_process(eps, theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15440 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 15004 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def pearson_correlation(x, y):\n",
    "    return np.mean((x - x.mean()) * (y - y.mean())) / (x.std() * y.std())\n",
    "\n",
    "def acf(x, lag=40):\n",
    "    \"\"\"\n",
    "    Determine auto correlation factors.\n",
    "    :param x: (array) Time series.\n",
    "    :param lag: (int) Number of lags.\n",
    "    \"\"\"\n",
    "    return np.array([1] + [pearson_correlation(x[:-i], x[i:]) for i in range(1, lag)])\n",
    "\n",
    "lag = 40\n",
    "# Create an ma(1) and an ma(2) process.\n",
    "ma_1 = ma_process(eps, [1])\n",
    "ma_2 = ma_process(eps, [0.2, -0.3, 0.8])\n",
    "\n",
    "\n",
    "\n",
    "def bartletts_formula(acf_array, n):\n",
    "    \"\"\"\n",
    "    Computes the Standard Error of an acf with Bartlet's formula\n",
    "    Read more at: https://en.wikipedia.org/wiki/Correlogram\n",
    "    :param acf_array: (array) Containing auto correlation factors\n",
    "    :param n: (int) Length of original time series sequence.\n",
    "    \"\"\"\n",
    "    # The first value has autocorrelation with it self. So that values is skipped\n",
    "    se = np.zeros(len(acf_array) - 1)\n",
    "    se[0] = 1 / np.sqrt(n)\n",
    "    se[1:] = np.sqrt((1 + 2 * np.cumsum(acf_array[1:-1]**2)) / n )\n",
    "    return se\n",
    "\n",
    "\n",
    "def plot_acf(x, alpha=0.05, lag=40, title=None):\n",
    "    \"\"\"\n",
    "    :param x: (array)\n",
    "    :param alpha: (flt) Statistical significance for confidence interval.\n",
    "    :parm lag: (int)\n",
    "    \"\"\"\n",
    "    acf_val = acf(x, lag)\n",
    "    plt.figure(figsize=(16, 4))\n",
    "    plt.vlines(np.arange(lag), 0, acf_val)\n",
    "    plt.scatter(np.arange(lag), acf_val, marker='o')\n",
    "    plt.xlabel('lag')\n",
    "    plt.ylabel('autocorrelation')\n",
    "    \n",
    "    # Determine confidence interval\n",
    "    ci = stats.norm.ppf(1 - alpha / 2.) * bartletts_formula(acf_val, len(x))\n",
    "    plt.fill_between(np.arange(1, ci.shape[0] + 1), -ci, ci, alpha=0.25)\n",
    "\n",
    "for array in [ma_1, ma_2]:\n",
    "    plot_acf(array)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 238816 characters]\n",
      "text/plain": [
       "<Figure size 1152x864 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def ar_process(eps, phi):\n",
    "    \"\"\"\n",
    "    Creates a AR process with a zero mean.\n",
    "    \"\"\"\n",
    "    # Reverse the order of phi and add a 1 for current eps_t\n",
    "    phi = np.r_[1, phi][::-1] \n",
    "    ar = eps.copy()\n",
    "    offset = len(phi)\n",
    "    for i in range(offset, ar.shape[0]):\n",
    "        ar[i - 1] = ar[i - offset: i] @ phi\n",
    "    return ar\n",
    "\n",
    "fig = plt.figure(figsize=(16, 4 * 3))\n",
    "a = 310\n",
    "for i in range(0, 11, 5):\n",
    "    a += 1\n",
    "    phi = np.random.normal(0, 0.1, size=i + 1)\n",
    "    plt.subplot(a)\n",
    "    plt.title(f'$\\\\phi$ = {phi.round(2)}')\n",
    "    plt.plot(ar_process(eps, phi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 16540 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 15688 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 17764 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_acf(ar_process(eps, [0.3, -0.3, 0.5]))\n",
    "plot_acf(ar_process(eps, [0.5, -0.1, 0.1]))\n",
    "plot_acf(ar_process(eps, [0.2, 0.5, 0.1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary least squares\n",
    "\n",
    "$$y = \\beta X + \\epsilon$$\n",
    "\n",
    "$$\\beta = (X^T X)^{-1} X^T y$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_squares(x, y):\n",
    "    return np.linalg.inv((x.T @ x)) @ (x.T @ y)\n",
    "\n",
    "class LinearModel:\n",
    "    def __init__(self, fit_intercept=True):\n",
    "        self.fit_intercept = fit_intercept\n",
    "        self.beta = None\n",
    "        self.intercept_ = None\n",
    "        self.coef_ = None\n",
    "    \n",
    "    def _prepare_features(self, x):\n",
    "        if self.fit_intercept:\n",
    "            x = np.hstack((np.ones((x.shape[0], 1)), x))\n",
    "        return x\n",
    "    \n",
    "    def fit(self, x, y):\n",
    "        x = self._prepare_features(x)\n",
    "        self.beta = least_squares(x, y)\n",
    "        if self.fit_intercept:\n",
    "            self.intercept_ = self.beta[0]\n",
    "            self.coef_ = self.beta[1:]\n",
    "        else:\n",
    "            self.coef_ = self.beta\n",
    "        \n",
    "    def predict(self, x):\n",
    "        x = self._prepare_features(x)\n",
    "        return x @ self.beta\n",
    "    \n",
    "    def fit_predict(self, x, y):\n",
    "        self.fit(x, y)\n",
    "        return self.predict(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15484 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 15644 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 15168 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "def pacf(x, lag=40):\n",
    "    \"\"\"\n",
    "    Partial auto correlation function.\n",
    "    \n",
    "    pacf results in:\n",
    "        [1, acf_lag_1, pacf_lag_2, pacf_lag_3]\n",
    "    :param x: (array)\n",
    "    :param lag: (int)\n",
    "    \"\"\"\n",
    "    y = []\n",
    "    \n",
    "    # Partial auto correlation needs intermediate terms.\n",
    "    # Therefore we start at index 3\n",
    "    for i in range(3, lag + 2):\n",
    "        backshifted = lag_view(x, i)[0]\n",
    "\n",
    "        xt = backshifted[:, 0]\n",
    "        feat = backshifted[:, 1:-1]\n",
    "        xt_hat = LinearModel(fit_intercept=False).fit_predict(feat, xt)\n",
    "\n",
    "        xt_k = backshifted[:, -1]\n",
    "        xt_k_hat = LinearModel(fit_intercept=False).fit_predict(feat, xt_k)\n",
    "\n",
    "        y.append(pearson_correlation(xt - xt_hat, xt_k - xt_k_hat))\n",
    "    return np.array([1, acf(x, 2)[1]] +  y)\n",
    "\n",
    "def plot_pacf(x, alpha=0.05, lag=40, title=None):\n",
    "    \"\"\"\n",
    "    :param x: (array)\n",
    "    :param alpha: (flt) Statistical significance for confidence interval.\n",
    "    :parm lag: (int)\n",
    "    \"\"\"\n",
    "    pacf_val = pacf(x, lag)\n",
    "    plt.figure(figsize=(16, 4))\n",
    "    plt.vlines(np.arange(lag + 1), 0, pacf_val)\n",
    "    plt.scatter(np.arange(lag + 1), pacf_val, marker='o')\n",
    "    plt.xlabel('lag')\n",
    "    plt.ylabel('autocorrelation')\n",
    "    \n",
    "    # Determine confidence interval\n",
    "    ci = stats.norm.ppf(1 - alpha / 2.) * bartletts_formula(pacf_val, len(x))\n",
    "    plt.fill_between(np.arange(1, ci.shape[0] + 1), -ci, ci, alpha=0.25)\n",
    "\n",
    "plot_pacf(ar_process(eps, [0.3, -0.3, 0.5]))\n",
    "plot_pacf(ar_process(eps, [0.5, -0.1, 0.1]))\n",
    "plot_pacf(ar_process(eps, [0.2, 0.5, 0.1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9, -1, -1, -1, -1, -1, -1, -1, -1, -1])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def difference(x, d=1):\n",
    "    if d == 0:\n",
    "        return x\n",
    "    else:\n",
    "        x = np.r_[x[0], np.diff(x)]\n",
    "        return difference(x, d - 1)\n",
    "    \n",
    "def undo_difference(x, d=1):\n",
    "    if d == 1:\n",
    "        return np.cumsum(x)\n",
    "    else:\n",
    "        x = np.cumsum(x)\n",
    "        return undo_difference(x, d - 1)\n",
    "    \n",
    "difference(np.arange(10)[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class AR(LinearModel):\n",
    "    def __init__(self, p, fit_intercept=True):\n",
    "        \"\"\"\n",
    "        An AR model.\n",
    "        :param p: (int)\n",
    "        :param fit_intercept: (bool)\n",
    "        \"\"\"\n",
    "        super().__init__(fit_intercept)\n",
    "        self.p = p\n",
    "        self.resid = None\n",
    "        self.variance = None\n",
    "        \n",
    "    def prepare_features(self, x):\n",
    "        features, x = lag_view(np.r_[np.zeros(self.p), x], self.p)\n",
    "        return features, x\n",
    "        \n",
    "    def fit_predict(self, x):\n",
    "        features, x = self.prepare_features(x)\n",
    "        super().fit(features, x)\n",
    "        return self.predict(x, prepared=(features, x))\n",
    "        \n",
    "    def predict(self, x, **kwargs):\n",
    "        \"\"\"\n",
    "        :param x: (array)\n",
    "        :kwargs:\n",
    "            prepared: (tpl) containing the features and x\n",
    "        \"\"\"\n",
    "        features, x = kwargs.get('prepared', self.prepare_features(x))\n",
    "        y = super().predict(features)\n",
    "        self.resid = x - y\n",
    "        self.variance = np.mean(self.resid**2)\n",
    "        return y\n",
    "    \n",
    "    @property\n",
    "    def bic(self):\n",
    "        n = len(self.resid)\n",
    "        return np.log(self.variance) + len(self.beta) * np.log(n) / n\n",
    "    \n",
    "    @property\n",
    "    def aic(self):\n",
    "        return np.log(self.variance) + 2 * len(self.beta) / len(self.resid)\n",
    "    \n",
    "# p = 1\n",
    "# m = AR(p)\n",
    "# pred = m.fit_predict(x)\n",
    "\n",
    "# print(m.coef_, m.intercept_)\n",
    "# plt.plot(x)\n",
    "# plt.plot(x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[19.43421933 -0.26098511 -0.23854199  1.10820997  0.20252782]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 81124 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class ARIMA(LinearModel):\n",
    "    def __init__(self, q, d, p):\n",
    "        \"\"\"\n",
    "        An ARIMA model.\n",
    "        :param q: (int) Order of the MA model.\n",
    "        :param p: (int) Order of the AR model.\n",
    "        :param d: (int) Number of times the data needs to be differenced.\n",
    "        \"\"\"\n",
    "        super().__init__(True)\n",
    "        self.p = p\n",
    "        self.d = d\n",
    "        self.q = q\n",
    "        self.ar = None\n",
    "        self.resid = None\n",
    "        \n",
    "    def prepare_features(self, x):\n",
    "        if self.d > 0:\n",
    "            x = difference(x, self.d)\n",
    "                    \n",
    "        ar_features = None\n",
    "        ma_features = None\n",
    "        \n",
    "        # Determine the features and the epsilon terms for the MA process\n",
    "        if self.q > 0:\n",
    "            if self.ar is None:\n",
    "                self.ar = ARIMA(0, 0, self.p)\n",
    "                self.ar.fit_predict(x)\n",
    "            eps = self.ar.resid\n",
    "            eps[0] = 0\n",
    "            \n",
    "            # prepend with zeros as there are no residuals_t-k in the first X_t\n",
    "            ma_features, _ = lag_view(np.r_[np.zeros(self.q), eps], self.q)\n",
    "            \n",
    "        # Determine the features for the AR process\n",
    "        if self.p > 0:\n",
    "            # prepend with zeros as there are no X_t-k in the first X_t\n",
    "            ar_features = lag_view(np.r_[np.zeros(self.p), x], self.p)[0]\n",
    "                                \n",
    "        if ar_features is not None and ma_features is not None:\n",
    "            n = min(len(ar_features), len(ma_features)) \n",
    "            ar_features = ar_features[:n]\n",
    "            ma_features = ma_features[:n]\n",
    "            features = np.hstack((ar_features, ma_features))\n",
    "        elif ma_features is not None: \n",
    "            n = len(ma_features)\n",
    "            features = ma_features[:n]\n",
    "        else:\n",
    "            n = len(ar_features)\n",
    "            features = ar_features[:n]\n",
    "        \n",
    "        return features, x[:n]\n",
    "    \n",
    "    def fit(self, x):\n",
    "        features, x = self.prepare_features(x)\n",
    "        super().fit(features, x)\n",
    "        return features\n",
    "            \n",
    "    def fit_predict(self, x): \n",
    "        \"\"\"\n",
    "        Fit and transform input\n",
    "        :param x: (array) with time series.\n",
    "        \"\"\"\n",
    "        features = self.fit(x)\n",
    "        return self.predict(x, prepared=(features))\n",
    "    \n",
    "    def predict(self, x, **kwargs):\n",
    "        \"\"\"\n",
    "        :param x: (array)\n",
    "        :kwargs:\n",
    "            prepared: (tpl) containing the features, eps and x\n",
    "        \"\"\"\n",
    "        features = kwargs.get('prepared', None)\n",
    "        if features is None:\n",
    "            features, x = self.prepare_features(x)\n",
    "        \n",
    "        y = super().predict(features)\n",
    "        self.resid = x - y\n",
    "\n",
    "        return self.return_output(y)\n",
    "    \n",
    "    def return_output(self, x):\n",
    "        if self.d > 0:\n",
    "            x = undo_difference(x, self.d) \n",
    "        return x\n",
    "    \n",
    "    def forecast(self, x, n):\n",
    "        \"\"\"\n",
    "        Forecast the time series.\n",
    "        \n",
    "        :param x: (array) Current time steps.\n",
    "        :param n: (int) Number of time steps in the future.\n",
    "        \"\"\"\n",
    "        features, x = self.prepare_features(x)\n",
    "        y = super().predict(features)\n",
    "        \n",
    "        # Append n time steps as zeros. Because the epsilon terms are unknown\n",
    "        y = np.r_[y, np.zeros(n)]\n",
    "        for i in range(n):\n",
    "            feat = np.r_[y[-(self.p + n) + i: -n + i], np.zeros(self.q)]\n",
    "            y[x.shape[0] + i] = super().predict(feat[None, :])\n",
    "        return self.return_output(y)\n",
    "    \n",
    "         \n",
    "q = 1\n",
    "d = 0\n",
    "p = 3\n",
    "\n",
    "\n",
    "m = ARIMA(q, d, p)\n",
    "pred = m.fit_predict(x)\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(x[a:b])\n",
    "ax.plot(pred[a:b])\n",
    "print(m.beta)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff87c89fe10>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 65752 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = sm.datasets.sunspots.load_pandas().data\n",
    "x = data['SUNACTIVITY'].values.squeeze()\n",
    "\n",
    "plt.figure(figsize=(16,4))\n",
    "plt.title('Sunspots')\n",
    "plt.plot(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20568 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 17356 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_acf(x)\n",
    "plot_pacf(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<__main__.ARIMA object at 0x7ff87a2ab438>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff87a2764e0>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 122888 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "q = 1\n",
    "d = 0\n",
    "p = 3\n",
    "\n",
    "m = ARIMA(q, d, p)\n",
    "pred = m.fit_predict(x)\n",
    "\n",
    "plt.figure(figsize=(16, 6))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(pred, label='forecast')\n",
    "ax.plot(x, label='y')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 125468 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = sm.tsa.ARIMA(x, (p, d, q)).fit()\n",
    "\n",
    "plt.figure(figsize=(16, 6))\n",
    "ax = plt.subplot(111)\n",
    "pred_sm = results.plot_predict(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff8785f83c8>]"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 78380 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pred = m.forecast(x, 40)\n",
    "plt.figure(figsize=(12,4))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(x[a:b])\n",
    "ax.plot(pred[a:b])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "::\n",
      "\n",
      "    Number of Observations - 309 (Annual 1700 - 2008)\n",
      "    Number of Variables - 1\n",
      "    Variable name definitions::\n",
      "\n",
      "        SUNACTIVITY - Number of sunspots for each year\n",
      "\n",
      "    The data file contains a 'YEAR' variable that is not returned by load.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sm.datasets.sunspots.NOTE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.sunspots.load_pandas().data\n",
    "data.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))\n",
    "del data['YEAR']\n",
    "x = data.values.squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 59164 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data.plot(figsize=(12, 4))\n",
    "x = data.values.squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moveing average model MA(q)\n",
    "\n",
    "$$ X_t = \\mu + \\epsilon_t + \\theta_1 \\epsilon_{t-1} ... + \\theta_q \\epsilon_{q} $$\n",
    "\n",
    "where $\\epsilon$ is white noise:\n",
    "\n",
    "$$\\epsilon \\sim N(0, \\alpha)$$\n",
    "\n",
    "$\\alpha$ is a hyperparameter\n",
    "\n",
    "\n",
    "\n",
    "# ARIMA model\n",
    "\n",
    "The mathematical definition of the ARIMA model is equal to that of the ARMA model.\n",
    "\n",
    "**ARMA(p, q)** _(AutoRegressive, Moving Average)_\n",
    "\n",
    "$$X_t = c + \\epsilon_t + \\sum_{i=1}^{p}{\\phi_i X_{t - i}} + \\sum_{i = 1}^q{\\theta_i \\epsilon_{t-i}} $$\n",
    "\n",
    "\n",
    "To convert this to an ARIMA model:\n",
    "\n",
    "**ARIMA(p, d, q)** _(AutoRegressive Integrated Moving Average)_\n",
    "\n",
    "$$X^*_t = c + \\epsilon_t + \\sum_{i=1}^{p}{\\phi_i X^*_{t - i}} + \\sum_{i = 1}^q{\\theta_i \\epsilon_{t-i}} $$\n",
    "\n",
    "Where $X^*$ is the integrated time series $X$ $d$ times\n",
    "\n",
    "The term integrated implies that in order to obtain a forecast of $X$, we have to sum up (or integrate over) the values of $X^*$ because $X^*$ are the differenced values of the original series $X$.If no differencing is involved, this model is called an Autoregressive Moving Average [ARMA(p,q)] with p and q retaining their original meaning and no d.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "def difference(x, d=1):\n",
    "    if d == 0:\n",
    "        return x\n",
    "    else:\n",
    "        x = np.concatenate((np.zeros(1), np.diff(x)))\n",
    "        return difference(x, d - 1)\n",
    "    \n",
    "def undo_difference(x, d=1):\n",
    "    if d == 1:\n",
    "        return np.cumsum(x)\n",
    "    else:\n",
    "        x = np.cumsum(x)\n",
    "        return undo_difference(x, d - 1)\n",
    "\n",
    "# print(difference(a, 2))\n",
    "\n",
    "# print(difference(a, 3))\n",
    "# print(undo_difference(difference(a, 3), 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.cumsum(difference(a, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "68.2529386667413 [-0.21150225 -0.26660066 -0.38434853 -0.36914932 -0.41741958 -0.36728083]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1d451a74a8>]"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 82472 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def init_eps(order, y, sigma=None):\n",
    "    sigma = y.std() / 3 if sigma is None else sigma\n",
    "    return np.random.normal(0, sigma, size=(y.shape[0]))\n",
    "\n",
    "def create_ar_features(x, order):\n",
    "    \"\"\"\n",
    "    Create the features for AR model.\n",
    "    \n",
    "    The features are [X_t-1, X_t-2 ... X_t-k] and the labels are X_t\n",
    "    So the features are shifted by one time step compared to the time step we want to predict.\n",
    "    \"\"\"\n",
    "    y = x.copy()\n",
    "    # Create features by shifting the window of `order` size by one step.\n",
    "    # This results in a 2D array [[t1, t2, t3], [t2, t3, t4], ... [t_k-2, t_k-1, t_k]]\n",
    "    x = np.array([y[-(i + order):][:order] for i in range(y.shape[0])])\n",
    "    \n",
    "    # Reverse the array as we started at the end and remove duplicates.\n",
    "    # Note that we truncate the features [order -1:] and the labels [order]\n",
    "    # This is the shifting of the features with one time step compared to the labels\n",
    "    x = np.stack(x)[::-1][order - 1: -1]\n",
    "    y = y[order:]\n",
    "    \n",
    "    return x, y\n",
    "\n",
    "def create_ma_features(eps, order, y=None):\n",
    "    eps = np.array([eps[-(i + order):][:order] for i in range(y.shape[0])])\n",
    "    eps = np.stack(eps)[::-1][order - 1:]\n",
    "    # eps_t is known when predicting x_t so the labels are in the same phase /lag\n",
    "    # as the features eps.\n",
    "    if y is not None:\n",
    "        y = y[order - 1:]    \n",
    "        return eps, y\n",
    "    else: \n",
    "        return eps\n",
    "\n",
    "\n",
    "# https://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model \n",
    "def ar_model(y, eps, order=2):\n",
    "    x, y = create_ra_features(y)    \n",
    "    \n",
    "    # Move eps_t to other side of equation\n",
    "    add = eps[:y.shape[0]]\n",
    "    y -= add\n",
    "    \n",
    "    m = linear_model.LinearRegression(fit_intercept=True)\n",
    "    m.fit(x, y)\n",
    "    \n",
    "    return m.predict(x) + add, m\n",
    "\n",
    "\n",
    "class ARIMA:\n",
    "    def __init__(self, p, d, q, sigma=None):\n",
    "        self.p = p\n",
    "        self.d = d\n",
    "        self.q = q\n",
    "        self.m = None\n",
    "        self.sigma = sigma\n",
    "        \n",
    "    def prepare_features(self, x):\n",
    "        if self.d > 0:\n",
    "            x = difference(x, self.d)\n",
    "        eps = init_eps(self.q, x)\n",
    "        \n",
    "        x_ = None\n",
    "        ma_features = ar_features = np.zeros((x.shape[0], 1))\n",
    "        if self.q > 0:\n",
    "            # Add features for `order` + 1. Later this will be eps_t\n",
    "            ma_features, x_ = create_ma_features(eps, self.q + 1, x)\n",
    "        if self.p > 0:\n",
    "            ar_features, x = create_ar_features(x, self.p)\n",
    "        x = x if x_ is None or len(x) < len(x_) else x_\n",
    "        del x_\n",
    "            \n",
    "        eps = eps[:x.shape[0]]\n",
    "        ma_features = ma_features[:x.shape[0]]\n",
    "        ar_features = ar_features[:x.shape[0]]\n",
    "        eps_t = ma_features[:, :1]\n",
    "        ma_features = ma_features[:, 1:]\n",
    "        return np.hstack((ma_features, ar_features)), eps_t.ravel(), x\n",
    "        \n",
    "    def fit_predict(self, x): \n",
    "        # move eps_t to left side of the equation\n",
    "        \n",
    "        features, eps_t, x = self.prepare_features(x)\n",
    "\n",
    "        x = x - eps_t\n",
    "        model = LinearModel()\n",
    "        model.fit(features, x)\n",
    "        self.m = model\n",
    "        return self.predict(x, prepared=(features, eps_t, x))\n",
    "    \n",
    "    def predict(self, x, **kwargs):\n",
    "        \"\"\"\n",
    "        :param x: (array)\n",
    "        :kwargs:\n",
    "            prepared: (tpl) containing the features, eps and x\n",
    "        \"\"\"\n",
    "        features, eps_t, x = kwargs.get('prepared', self.prepare_features(x))\n",
    "\n",
    "        x = self.m.predict(features) + eps_t\n",
    "        if self.d > 0:\n",
    "            x += undo_difference(x, self.d)\n",
    "        return x\n",
    "    \n",
    "    def generate(self, n, ma_w=None, ar_w=None):\n",
    "        \n",
    "        t = np.arange(n)\n",
    "        if self.q:\n",
    "            eps = np.random.normal(0, 1, size=n)\n",
    "            feat = create_ma_features(eps)\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "q = 5\n",
    "d = 0\n",
    "p = 1\n",
    "\n",
    "m = ARIMA(p, d, q, sigma=30)\n",
    "pred = m.fit_predict(x)\n",
    "print(m.m.intercept_, m.m.coef_)\n",
    "# print(m.RA.intercept_)\n",
    "# m.RA.coef_\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(data)\n",
    "plt.plot(data.index[data.shape[0] - pred.shape[0]:], pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/statsmodels/tsa/base/tsa_model.py:171: ValueWarning: No frequency information was provided, so inferred frequency A-DEC will be used.\n",
      "  % freq, ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f76f1fe3a58>]"
      ]
     },
     "execution_count": 237,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 80860 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = 5\n",
    "q = 2\n",
    "d = 0\n",
    "results = sm.tsa.ARIMA(data.copy(), (p, d, q)).fit()\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(data)\n",
    "\n",
    "pred = results.model.predict(results.params, dynamic=False)\n",
    "\n",
    "plt.plot(data.index[:pred.shape[0]], pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f76f1fecac8>]"
      ]
     },
     "execution_count": 238,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 35524 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(results.resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Autocorrelation\n",
    "\n",
    "Let's look how autocorrelation is defined. First we make a function that computes the Pearson correlation between two sequences that are assumed to come from a Normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33916 characters]\n",
      "text/plain": [
       "<Figure size 864x576 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = fig.add_subplot(211)\n",
    "fig = sm.graphics.tsa.plot_acf(data.values.squeeze(), lags=40, ax=ax)\n",
    "ax = fig.add_subplot(212)\n",
    "fig = sm.graphics.tsa.plot_pacf(data.values.squeeze(), lags=40, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9996781213197145"
      ]
     },
     "execution_count": 240,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pearson_correlation(x, y):\n",
    "    return np.mean((x - x.mean()) * (y - y.mean())) / (x.std() * y.std())\n",
    "    \n",
    "pearson_correlation(x, x + np.random.randn(x.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The autocorrelation is the correlation between the sequence with the current last value truncated at the end `x[:-i]` and the sequence with a certain amount of lag `x[:i]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 8360 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot(x, lag=40):\n",
    "    plt.figure(figsize=(12, 4))\n",
    "    corr = [pearson_correlation(x[:-i], x[i:]) for i in range(1, lag)]\n",
    "    plt.vlines(np.arange(1, lag), 0, corr)\n",
    "    plt.show()\n",
    "plot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partial autocorrelation\n",
    "\n",
    "The partial correlation between $BX$ and $B^kX$  where $B$ is the backshift operator.\n",
    "\n",
    "$$BX_t = X_{t-1} $$\n",
    "\n",
    "can be determined as follows.\n",
    "\n",
    "Let $\\hat{BX}$ and $\\hat{B^kX}$ be determined by a Linear Model conditioned on $ B^2X ... B^{k -1} $ parameterized by $\\beta$ and $\\phi$.\n",
    "\n",
    "$$\\hat{BX} \\sim \\beta_2 B^2X + \\beta_3 B^3X ... \\beta_{k-1} B^{k -1}X $$ \n",
    "\n",
    "$$\\hat{B^kX} \\sim \\phi_2 B^2X + \\phi_3 B^3X ... \\phi_{k-1} B^{k -1}X $$ \n",
    "\n",
    "The partial correlation is then defined by the Pearsons coefficient of the residuals of both predicted values $\\hat{BX}$ and $\\hat{B^kX}$.\n",
    "\n",
    "$$ pcf = corr(BX - \\hat{BX} , B^kX - \\hat{B^kX} ) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5., 11., 16., 23.],\n",
       "       [11., 16., 23., 36.],\n",
       "       [16., 23., 36., 58.]])"
      ]
     },
     "execution_count": 242,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "backshifted = create_ar_features(x, 4)[0]\n",
    "backshifted[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Correlation of $X_{t-1}$, $X_{t-2} ... X_{t-n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcaf(x, lag=40):\n",
    "    y = []\n",
    "    for i in range(3, lag):\n",
    "        backshifted = create_ar_features(x, i)[0]\n",
    "\n",
    "        bx = backshifted[:, 0]\n",
    "        feat = backshifted[:, 1:-1]\n",
    "        bx_hat = LinearModel(fit_intercept=True).fit_predict(feat, bx)\n",
    "\n",
    "        bkx = backshifted[:, -1]\n",
    "        bkx_hat = LinearModel(fit_intercept=True).fit_predict(feat, bkx)\n",
    "\n",
    "        y.append(pearson_correlation(bx - bx_hat, bkx - bkx_hat))\n",
    "    return y\n",
    "\n",
    "y = pcaf(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 9868 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_autocorrelation(x, lag=40, biased=False):\n",
    "    y = []\n",
    "    for i in range(3, lag):\n",
    "        backshifted = create_ar_features(x, i)[0]\n",
    "\n",
    "        bx = backshifted[:, 0]\n",
    "        feat = backshifted[:, 1:-1]\n",
    "        bx_hat = LinearModel(fit_intercept=False).fit_predict(feat, bx)\n",
    "\n",
    "        bkx = backshifted[:, -1]\n",
    "        bkx_hat = LinearModel(fit_intercept=False).fit_predict(feat, bkx)\n",
    "\n",
    "        y.append(pearson_correlation(bx - bx_hat, bkx - bkx_hat))\n",
    "\n",
    "    plt.figure(figsize=(12, 4))\n",
    "    plt.ylim(-1, 1)\n",
    "    plt.vlines(np.arange(0, len(y)) + 2, 0, y)\n",
    "plot_autocorrelation(x, 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.   ,  0.823, -0.69 , -0.141,  0.056,  0.005,  0.182,  0.22 ,\n",
       "        0.228,  0.263, -0.016, -0.008, -0.015, -0.001,  0.062, -0.083,\n",
       "       -0.079, -0.162, -0.091,  0.045,  0.002,  0.113, -0.001, -0.099,\n",
       "       -0.057,  0.023, -0.066,  0.097,  0.111, -0.126,  0.025, -0.048,\n",
       "       -0.011,  0.005,  0.009,  0.02 , -0.048,  0.069, -0.041,  0.045,\n",
       "        0.048])"
      ]
     },
     "execution_count": 252,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16844 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12, 4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_ylim(-1, 1)\n",
    "sm.graphics.tsa.plot_pacf(x, lags=40, ax=ax, method='ols')\n",
    "\n",
    "ax.vlines(np.arange(0, len(y)) + 2, 0, y, color='r')\n",
    "\n",
    "sm.tsa.pacf(data).round(3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "|   :) | AR(p)                  | MA(q)                  | ARMA(p, q) |\n",
    "|------|----------------------- |------------------------|------------|\n",
    "| ACF  | Tails off              | Cuts off after lag $q$ | Tails off  |\n",
    "| PACF | Cuts off after lag $p$ | Tails off              | Tails off  |\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
