{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python [conda env:tensorflow2]",
      "language": "python",
      "name": "conda-env-tensorflow2-py"
    },
    "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.10"
    },
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Sahel-Eskandar/Project-Notebooks/blob/main/6_Anomaly_Detection_VAR.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9QXuZifmYRWN"
      },
      "source": [
        "#!pip install -U statsmodels"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZB_30mP2Wmpy"
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scipy.stats as stats\n",
        "import matplotlib.pyplot as plt\n",
        "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
        "from statsmodels.tsa.vector_ar.var_model import VAR\n",
        "\n",
        "import warnings\n",
        "warnings.simplefilter('ignore')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "D1NeoT3wWy1c",
        "outputId": "e97918c1-a868-4dbd-ed20-d7951ebcd2a4"
      },
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 212
        },
        "id": "p6cxykziWmp1",
        "outputId": "1ae72a83-66cb-4d90-cf85-c108e550f74d"
      },
      "source": [
        "### READ DATA ###\n",
        "\n",
        "train_hours = 80*7*24  # weeks x hours x days\n",
        "test_hours = 15*7*24  # weeks x hours x days\n",
        "\n",
        "df = pd.read_csv('/content/drive/MyDrive/VAR/bike-and-ped-counter.csv.zip', nrows=train_hours+test_hours, parse_dates=['Date'])\n",
        "\n",
        "print(df.shape)\n",
        "df.head()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(15960, 6)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Date</th>\n",
              "      <th>BGT North of NE 70th Total</th>\n",
              "      <th>Ped South</th>\n",
              "      <th>Ped North</th>\n",
              "      <th>Bike North</th>\n",
              "      <th>Bike South</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>2014-01-01 00:00:00</td>\n",
              "      <td>15.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.0</td>\n",
              "      <td>2.0</td>\n",
              "      <td>11.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>2014-01-01 01:00:00</td>\n",
              "      <td>9.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>7.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>2014-01-01 02:00:00</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>9.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>2014-01-01 03:00:00</td>\n",
              "      <td>19.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>19.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>2014-01-01 04:00:00</td>\n",
              "      <td>19.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>19.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                 Date  BGT North of NE 70th Total  ...  Bike North  Bike South\n",
              "0 2014-01-01 00:00:00                        15.0  ...         2.0        11.0\n",
              "1 2014-01-01 01:00:00                         9.0  ...         1.0         7.0\n",
              "2 2014-01-01 02:00:00                         9.0  ...         0.0         9.0\n",
              "3 2014-01-01 03:00:00                        19.0  ...         0.0        19.0\n",
              "4 2014-01-01 04:00:00                        19.0  ...         0.0        19.0\n",
              "\n",
              "[5 rows x 6 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tSxkOyJMWmp2"
      },
      "source": [
        "### FILL NAN ###\n",
        "\n",
        "df['Ped South'] = df['Ped South'].groupby(df.Date.dt.hour).transform(lambda x: x.fillna(x.median()))\n",
        "df['Ped North'] = df['Ped North'].groupby(df.Date.dt.hour).transform(lambda x: x.fillna(x.median()))\n",
        "df['Bike South'] = df['Bike South'].groupby(df.Date.dt.hour).transform(lambda x: x.fillna(x.median()))\n",
        "df['Bike North'] = df['Bike North'].groupby(df.Date.dt.hour).transform(lambda x: x.fillna(x.median()))\n",
        "df['BGT North of NE 70th Total'] = df['Ped South'] + df['Ped North'] + df['Bike South'] + df['Bike North']\n",
        "\n",
        "df['Date'] = pd.to_datetime(df['Date'].dt.date)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vi9Y6pfpWmp2"
      },
      "source": [
        "### DAILY AGGREGATION ###\n",
        "\n",
        "df_day = pd.DataFrame()\n",
        "df_day['Ped South'] = df.groupby(df.Date)['Ped South'].sum()\n",
        "df_day['Ped North'] = df.groupby(df.Date)['Ped North'].sum()\n",
        "df_day['Bike South'] = df.groupby(df.Date)['Bike South'].sum()\n",
        "df_day['Bike North'] = df.groupby(df.Date)['Bike North'].sum()\n",
        "df_day['Total'] = df.groupby(df.Date)['BGT North of NE 70th Total'].sum()\n",
        "\n",
        "df_day.index = pd.DatetimeIndex(df_day.index.values, freq=df_day.index.inferred_freq)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 403
        },
        "id": "GPxGU3MHWmp3",
        "outputId": "990199fb-0e55-468f-81b7-c3213d050833"
      },
      "source": [
        "### PLOT DAILY DATA ###\n",
        "\n",
        "df_day.plot(figsize=(16,6))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f86ed807da0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 207736 characters]\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FSE0_8ISWmp3",
        "outputId": "0fe3582d-81df-4bde-e8ff-b243a8f11114"
      },
      "source": [
        "### TRAIN SPLIT ###\n",
        "\n",
        "train = df_day[:(train_hours//24)].copy()\n",
        "\n",
        "train.shape"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(560, 5)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 409
        },
        "id": "AA-IlKsUWmp4",
        "outputId": "f33eae04-7ae0-41c7-c518-9b9687bcaede"
      },
      "source": [
        "### TOTAL COUNT TRAIN AUTOCORR ###\n",
        "\n",
        "plt.figure(figsize=(16,6))\n",
        "pd.plotting.autocorrelation_plot(train['Total'])"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f86ec3f76a0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 50304 characters]\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OqPhUWPZWmp4"
      },
      "source": [
        "### MONTHLY TRAIN MEAN ###\n",
        "\n",
        "month_mean_train = train.groupby(train.index.month).mean()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8FtL2M7IWmp4"
      },
      "source": [
        "### REMOVE LONG TERM SEASONALITY FROM TRAIN ###\n",
        "\n",
        "train['Ped South'] = train.apply(lambda x: x['Ped South'] - month_mean_train['Ped South'][x.name.month], axis=1)\n",
        "train['Ped North'] = train.apply(lambda x: x['Ped North'] - month_mean_train['Ped North'][x.name.month], axis=1)\n",
        "train['Bike South'] = train.apply(lambda x: x['Bike South'] - month_mean_train['Bike South'][x.name.month], axis=1)\n",
        "train['Bike North'] = train.apply(lambda x: x['Bike North'] - month_mean_train['Bike North'][x.name.month], axis=1)\n",
        "train['Total'] = train.apply(lambda x: x['Total'] - month_mean_train['Total'][x.name.month], axis=1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 409
        },
        "id": "CKV2QFV-Wmp5",
        "outputId": "16118c9d-4a71-4aa3-8a36-6eb31585952a"
      },
      "source": [
        "### TOTAL COUNT TRAIN AUTOCORR ###\n",
        "\n",
        "plt.figure(figsize=(16,6))\n",
        "pd.plotting.autocorrelation_plot(train['Total'])"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f86ebf90160>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 12
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 64924 characters]\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "14M1VBUqWmp5"
      },
      "source": [
        "# UNIVARIATE ANOMALY DETECTION"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bpgybS5uWmp6"
      },
      "source": [
        "### TRAIN TEST SPLIT IN UNIVARIATE CASE ###\n",
        "\n",
        "train_uni = train['Total'].copy()\n",
        "test_uni = df_day['Total'][(train_hours//24):].copy()\n",
        "test_uni = test_uni - test_uni.index.month.map(month_mean_train['Total'])\n",
        "\n",
        "train.drop('Total', inplace=True, axis=1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vcEJ0LfqWmp6",
        "outputId": "fb08bb5f-0362-40fe-b9eb-7612cc5639d8"
      },
      "source": [
        "### FIND BEST UNIVARIATE MODEL ###\n",
        "\n",
        "AIC = {}\n",
        "best_aic, best_order = np.inf, 0\n",
        "\n",
        "for p in range(6,9):\n",
        "    for q in range(0,10):\n",
        "\n",
        "        mod = SARIMAX(train_uni, order=(p,0,q), enforce_invertibility=False)\n",
        "        try:\n",
        "            res = mod.fit(disp=False)\n",
        "            AIC[(p,0,q)] = res.aic\n",
        "        except:\n",
        "            AIC[(p,0,q)] = np.inf\n",
        "            \n",
        "        if AIC[(p,0,q)] < best_aic:\n",
        "            best_aic = AIC[(p,0,q)]\n",
        "            best_order = (p,0,q)\n",
        "            \n",
        "print('BEST ORDER', best_order, 'BEST AIC:', best_aic)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "BEST ORDER (6, 0, 5) BEST AIC: 8525.504094721367\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pn4HJIkvWmp7",
        "outputId": "75e26d38-2c75-4ca6-a78b-3e37cbaa0f10"
      },
      "source": [
        "### TRAIN BEST UNIVARIATE MODEL ###\n",
        "\n",
        "mod = SARIMAX(train_uni, order=best_order, enforce_invertibility=False)\n",
        "res = mod.fit(disp=False)\n",
        "res.aic"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "8525.504094721367"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 621
        },
        "id": "F_OxObntWmp8",
        "outputId": "b266e8c8-e531-42cb-c28f-1eb206cdfbd3"
      },
      "source": [
        "### PERFORMANCE OPTINAL MODEL ###\n",
        "\n",
        "res.plot_diagnostics(figsize=(18,10))\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 113748 characters]\n",
            "text/plain": [
              "<Figure size 1296x720 with 4 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TmEZWIZtWmp8"
      },
      "source": [
        "### GET TRAIN FITTED VALUES ###\n",
        "\n",
        "predict = res.get_prediction()\n",
        "predicted_mean = predict.predicted_mean + predict.predicted_mean.index.month.map(month_mean_train['Total'])\n",
        "train_uni = train_uni + train_uni.index.month.map(month_mean_train['Total'])\n",
        "predict_ci = predict.conf_int(alpha=0.1)\n",
        "predict_ci['lower Total'] = predict_ci.apply(lambda x: x['lower Total'] + month_mean_train['Total'][x.name.month], axis=1)\n",
        "predict_ci['upper Total'] = predict_ci.apply(lambda x: x['upper Total'] + month_mean_train['Total'][x.name.month], axis=1)\n",
        "\n",
        "residuals_mean = res.resid.mean()\n",
        "residuals_std = res.resid.std()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_esjGh42Wmp9",
        "outputId": "04b3bce1-4e42-437b-9207-92c08edd16f3"
      },
      "source": [
        "### PLOT TRAIN FITTED VALUES ###\n",
        "\n",
        "plt.figure(figsize=(15,6))\n",
        "\n",
        "plt.plot(train_uni, 'k.')\n",
        "plt.plot(predicted_mean, linestyle='--', linewidth=2, color='blue')\n",
        "plt.fill_between(predict_ci.index, predict_ci['lower Total'], predict_ci['upper Total'], alpha=0.6)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PolyCollection at 0x7f86e66f59b0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 18
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 177704 characters]\n",
            "text/plain": [
              "<Figure size 1080x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y3iwZ7ktWmp9"
      },
      "source": [
        "### ITERATIVE PREDICTIONS ON TEST DATA ###\n",
        "\n",
        "point_forecast = res.get_prediction(end=mod.nobs)\n",
        "point_ci = point_forecast.conf_int(alpha=0.1)\n",
        "\n",
        "mean_pred = {point_forecast.predicted_mean.index[-1]: point_forecast.predicted_mean[-1]}\n",
        "upper_pred = {point_ci.index[-1]: point_ci['upper Total'][-1]}\n",
        "lower_pred = {point_ci.index[-1]: point_ci['lower Total'][-1]}\n",
        "\n",
        "for t,row in test_uni[:-1].iteritems():\n",
        "    \n",
        "    row = pd.Series(row, index=[t])\n",
        "    res = res.extend(row)\n",
        "    point_forecast = res.get_prediction(1)\n",
        "    point_ci = point_forecast.conf_int(alpha=0.1)\n",
        "    \n",
        "    mean_pred[point_forecast.predicted_mean.index[0]] = point_forecast.predicted_mean.values[0]\n",
        "    upper_pred[point_ci.index[0]] = point_ci['upper Total'][0]\n",
        "    lower_pred[point_ci.index[0]] = point_ci['lower Total'][0]\n",
        "    \n",
        "mean_pred = pd.Series(mean_pred)\n",
        "upper_pred = pd.Series(upper_pred)\n",
        "lower_pred = pd.Series(lower_pred)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "mGcnBdA7Wmp9",
        "outputId": "e0d97038-457e-43a3-c13c-39326cd1d435"
      },
      "source": [
        "### TEST RESIDUALS ANALYSIS ###\n",
        "\n",
        "alpha = 0.01\n",
        "upper = stats.norm.ppf(1 - alpha/2)\n",
        "lower = stats.norm.ppf(alpha/2)\n",
        "\n",
        "residuals_test = test_uni - residuals_mean\n",
        "residuals_test = (residuals_test - residuals_mean) / residuals_std\n",
        "\n",
        "plt.figure(figsize=(15,6))\n",
        "plt.plot(residuals_test)\n",
        "plt.ylabel('resid')\n",
        "plt.axhline(upper, c='red', linestyle='--')\n",
        "plt.axhline(lower, c='red', linestyle='--')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f86e66ac5c0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 20
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 61564 characters]\n",
            "text/plain": [
              "<Figure size 1080x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iwIYNpfnWmp-"
      },
      "source": [
        "### REVERSE SCALING TEST ITERATIVE PREDICTIONS ###\n",
        "\n",
        "mean_pred = mean_pred + mean_pred.index.month.map(month_mean_train['Total'])\n",
        "upper_pred = upper_pred + upper_pred.index.month.map(month_mean_train['Total'])\n",
        "lower_pred = lower_pred + lower_pred.index.month.map(month_mean_train['Total'])\n",
        "test_uni = test_uni + test_uni.index.month.map(month_mean_train['Total'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wby6PPN4Wmp-",
        "outputId": "31cc1d24-7e4d-4485-9d46-6d0ecbc641b1"
      },
      "source": [
        "### PLOT TEST ITERATIVE PREDICTIONS ###\n",
        "\n",
        "plt.figure(figsize=(15,6))\n",
        "\n",
        "plt.plot(test_uni, 'k.')\n",
        "plt.plot(mean_pred, linestyle='--', linewidth=2, color='blue')\n",
        "plt.fill_between(mean_pred.index, lower_pred, upper_pred, alpha=0.6)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PolyCollection at 0x7f86e6605ef0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 22
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 85408 characters]\n",
            "text/plain": [
              "<Figure size 1080x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yVZgjDzIWmp_"
      },
      "source": [
        "# MULTIVARIATE ANOMALY DETECTION"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "olBpmRCaWmp_",
        "outputId": "c6029afc-56ee-4c1b-f0e2-3eae90d8cc2d"
      },
      "source": [
        "### FIND BEST MULTIVARIATE MODEL ###\n",
        "\n",
        "AIC = {}\n",
        "best_aic, best_order = np.inf, 0\n",
        "\n",
        "for i in range(1,50):\n",
        "    model = VAR(endog=train)\n",
        "    var_result = model.fit(maxlags=i)\n",
        "    AIC[i] = var_result.aic\n",
        "    \n",
        "    if AIC[i] < best_aic:\n",
        "        best_aic = AIC[i]\n",
        "        best_order = i\n",
        "        \n",
        "print('BEST ORDER', best_order, 'BEST AIC:', best_aic)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "BEST ORDER 7 BEST AIC: 34.896281239608705\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CzGwGabyWmp_",
        "outputId": "15082027-987b-4f84-aa16-a0696f75dc8b"
      },
      "source": [
        "### TRAIN BEST MULTIVARIATE MODEL ###\n",
        "\n",
        "var = VAR(endog=train)\n",
        "var_result = var.fit(maxlags=best_order)\n",
        "\n",
        "var_result.aic"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "34.896281239608705"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 24
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7bAbCBGIWmp_"
      },
      "source": [
        "### COMPUTE TRAIN T2 METRIC ###\n",
        "\n",
        "residuals_mean = var_result.resid.values.mean(axis=0)\n",
        "residuals_std = var_result.resid.values.std(axis=0)\n",
        "\n",
        "residuals = (var_result.resid.values - residuals_mean) / residuals_std\n",
        "cov_residuals = np.linalg.inv(np.cov(residuals.T))\n",
        "\n",
        "T = np.diag((residuals).dot(cov_residuals).dot(residuals.T))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SozCFOB6WmqA",
        "outputId": "11ce9dd4-e348-4934-a7fe-8049f0f3b42f"
      },
      "source": [
        "### COMPUTE UCL ###\n",
        "\n",
        "m = var_result.nobs\n",
        "p = var_result.resid.shape[-1]\n",
        "alpha = 0.01\n",
        "\n",
        "UCL = stats.f.ppf(1-alpha, dfn=p, dfd=m-p) *(p*(m+1)*(m-1)/(m*m-m*p))\n",
        "UCL"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "13.511575593053317"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 26
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BSKfpfWUWmqA",
        "outputId": "6ca8a72f-f4b9-4665-9461-52b2b10d0ffc"
      },
      "source": [
        "### PLOT TRAIN T2 METRIC ###\n",
        "\n",
        "plt.figure(figsize=(16,6))\n",
        "plt.plot(train.index[best_order:], T)\n",
        "plt.ylabel('T-squared')\n",
        "plt.axhline(UCL, c='red', linestyle='--')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f86e656a208>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 27
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 39196 characters]\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-q3ltyWWWmqA",
        "outputId": "6f44107f-1591-4ad2-f9d8-797333f08346"
      },
      "source": [
        "### TEST SPLIT IN MULTIVARIATE CASE ###\n",
        "\n",
        "test = df_day[(train_hours//24-best_order):].copy()\n",
        "test.drop('Total', inplace=True, axis=1)\n",
        "\n",
        "test.shape"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(112, 4)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 28
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CpaAZW5wWmqB"
      },
      "source": [
        "### REMOVE LONG TERM SEASONALITY FROM TEST ###\n",
        "\n",
        "test['Ped South'] = test.apply(lambda x: x['Ped South'] - month_mean_train['Ped South'][x.name.month], axis=1)\n",
        "test['Ped North'] = test.apply(lambda x: x['Ped North'] - month_mean_train['Ped North'][x.name.month], axis=1)\n",
        "test['Bike South'] = test.apply(lambda x: x['Bike South'] - month_mean_train['Bike South'][x.name.month], axis=1)\n",
        "test['Bike North'] = test.apply(lambda x: x['Bike North'] - month_mean_train['Bike North'][x.name.month], axis=1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Yh_AznfkWmqB",
        "outputId": "e3c34ed3-b173-4803-d6e4-1952fb40a007"
      },
      "source": [
        "### ITERATIVE PREDICTIONS ON TEST DATA ###\n",
        "\n",
        "pred = []\n",
        "\n",
        "for i in range(best_order, len(test)):\n",
        "    \n",
        "    pred.append(var_result.forecast(test.iloc[i-best_order:i].values, steps=1))\n",
        "    \n",
        "pred = np.vstack(pred)\n",
        "pred.shape"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(105, 4)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 30
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qFqHDt3UWmqB"
      },
      "source": [
        "### COMPUTE TEST T2 METRIC ###\n",
        "\n",
        "residuals_test = test.iloc[best_order:].values - pred\n",
        "residuals_test = (residuals_test - residuals_mean) / residuals_std\n",
        "\n",
        "T_test = np.diag((residuals_test).dot(cov_residuals).dot(residuals_test.T))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "bGNimeG8WmqB",
        "outputId": "1c177c4e-50e5-45d2-a977-52ded87a6d1c"
      },
      "source": [
        "### PLOT TEST T2 METRIC ###\n",
        "\n",
        "plt.figure(figsize=(16,6))\n",
        "plt.plot(test.iloc[best_order:].index, T_test)\n",
        "plt.ylabel('T-squared')\n",
        "plt.axhline(UCL, c='red', linestyle='--')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f86e6544be0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 32
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "[encoded data, 38612 characters]\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}