10_VAR_VARMA_VARMAX_starter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "94b82def",
   "metadata": {},
   "source": [
    "# VAR, VARMA and VARMAX "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ef6b340",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.stattools import grangercausalitytests, adfuller\n",
    "from statsmodels.stats.diagnostic import acorr_ljungbox\n",
    "from statsmodels.tsa.statespace.varmax import VARMAX\n",
    "from tqdm import tqdm_notebook\n",
    "from itertools import product\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33cdb7eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (9,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e03ee61",
   "metadata": {},
   "source": [
    "Source: Statistics Canada. Table 32-10-0077-01  [Farm product prices, crops and livestock](https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=3210007701)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7333fc7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5f9214",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=1, dpi=300, figsize=(12,9))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    data = df[df.columns[i+1]]\n",
    "    \n",
    "    ax.plot(data, color='black', linewidth=1)\n",
    "    ax.set_title(df.columns[i+1])\n",
    "    ax.xaxis.set_ticks_position('none')\n",
    "    ax.yaxis.set_ticks_position('none')\n",
    "    ax.spines['top'].set_alpha(0)\n",
    "    ax.tick_params(labelsize=6)\n",
    "\n",
    "plt.setp(axes, xticks=np.arange(0, 369, 24), xticklabels=np.arange(1992, 2023, 2))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d201ff4e",
   "metadata": {},
   "source": [
    "## VAR(p) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cb04a55",
   "metadata": {},
   "source": [
    "### Check for stationarity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f0cd46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the ADF test for the \"cows\" series\n",
    "\n",
    "print('cows')\n",
    "print(f'ADF Statistic: {ad_fuller_result_1[0]}')\n",
    "print(f'p-value: {ad_fuller_result_1[1]}')\n",
    "\n",
    "print('\\n---------------------\\n')\n",
    "\n",
    "# Run the ADF test for the \"calves\" series\n",
    "\n",
    "print('calves')\n",
    "print(f'ADF Statistic: {ad_fuller_result_2[0]}')\n",
    "print(f'p-value: {ad_fuller_result_2[1]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c372c931",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the ADF test for the differenced \"cows\" series\n",
    "\n",
    "print('cows')\n",
    "print(f'ADF Statistic: {ad_fuller_result_1[0]}')\n",
    "print(f'p-value: {ad_fuller_result_1[1]}')\n",
    "\n",
    "print('\\n---------------------\\n')\n",
    "\n",
    "# Run the ADF test for the differences \"calves\" series\n",
    "\n",
    "print('calves')\n",
    "print(f'ADF Statistic: {ad_fuller_result_2[0]}')\n",
    "print(f'p-value: {ad_fuller_result_2[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab3573e",
   "metadata": {},
   "source": [
    "### Model selection "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f353e40",
   "metadata": {},
   "outputs": [],
   "source": [
    "def VAR_gridsearch(endog, min_p, max_p):\n",
    "    \n",
    "    # Determine the range of possible values for p\n",
    "    \n",
    "    # Initialize an empty list to store the order and AIC\n",
    "    \n",
    "    print(f'Fitting {all_p} unique models')\n",
    "    \n",
    "    for p in tqdm_notebook(all_p):\n",
    "        try:\n",
    "            # Fit a VAR(p) model\n",
    "        except:\n",
    "            continue\n",
    "        \n",
    "        # Append the order and AIC to the list\n",
    "\n",
    "    # Store the results in a DataFrame and rename the columns\n",
    "        \n",
    "    # Sort the DataFrame in ascending order of AIC\n",
    "    \n",
    "    return result_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e86ce30a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the endogenous variables\n",
    "\n",
    "# Get the differenced endogenous variables\n",
    "\n",
    "# Train/test split (last 48 times steps are for the test set)\n",
    "\n",
    "# Run the grid search\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f5aa6e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('calves Granger-causes cows?\\n')\n",
    "print('------------------')\n",
    "# Run the Granger causality test calves -> cows\n",
    "\n",
    "print('\\ncows Granger-causes calves?\\n')\n",
    "print('------------------')\n",
    "# Run the Granger causality test cows -> calves\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "565f1e53",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e9d3267b",
   "metadata": {},
   "source": [
    "### Residuals analysis\n",
    "#### cows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c43a501",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cows\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccf72b4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ljung_box_test(residuals, is_seasonal, period):\n",
    "    \n",
    "    if is_seasonal:\n",
    "        lb_df = acorr_ljungbox(residuals, period=period)\n",
    "    else:\n",
    "        max_lag = min([10, len(residuals)/5])\n",
    "        \n",
    "        lb_df = acorr_ljungbox(residuals, np.arange(1, max_lag+1, 1))\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(lb_df['lb_pvalue'], 'b-', label='p-values')\n",
    "    ax.hlines(y=0.05, xmin=1, xmax=len(lb_df), color='black')\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    if all(pvalue > 0.05 for pvalue in lb_df['lb_pvalue']):\n",
    "        print('All values are above 0.05. We fail to reject the null hypothesis. The residuals are uncorrelated')\n",
    "    else:\n",
    "        print('At least one p-value is smaller than 0.05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b08d7e91",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "59305e43",
   "metadata": {},
   "source": [
    "#### calves "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca0094a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calves\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8bf26fd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "700c0e56",
   "metadata": {},
   "source": [
    "### Forecasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9debb1b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rolling_predictions(df, last_train_value, train_len, horizon, window, method):\n",
    "    \n",
    "    total_len = train_len + horizon\n",
    "    \n",
    "    if method == 'VAR':\n",
    "\n",
    "        # Get the predictions from the VAR model\n",
    "        \n",
    "        return cows_pred_VAR[:horizon], calves_pred_VAR[:horizon]\n",
    "    \n",
    "    elif method == 'last':\n",
    "        cows_pred_last = []\n",
    "        calves_pred_last = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            \n",
    "            cows_last = df[:i].iloc[-1]['cows']\n",
    "            calves_last = df[:i].iloc[-1]['calves']\n",
    "            \n",
    "            cows_pred_last.extend(cows_last for _ in range(window))\n",
    "            calves_pred_last.extend(calves_last for _ in range(window))\n",
    "        \n",
    "        cows_pred_last = np.insert(cows_pred_last, 0, last_train_value['cows'])\n",
    "        cows_pred_last = cows_pred_last.cumsum()\n",
    "        \n",
    "        calves_pred_last = np.insert(calves_pred_last, 0, last_train_value['calves'])\n",
    "        calves_pred_last = calves_pred_last.cumsum()\n",
    "            \n",
    "        return cows_pred_last[:horizon], calves_pred_last[:horizon]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23cada9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_df = df[['cows', 'calves']][-48:]\n",
    "\n",
    "TRAIN_LEN = len(train)\n",
    "HORIZON = len(test)\n",
    "LAST_TRAIN_VALUE = endog.iloc[321]\n",
    "\n",
    "windows = [1, 4, 6, 12]\n",
    "\n",
    "for window in windows:\n",
    "    cows_pred_last, calves_pred_last = rolling_predictions(endog_diff, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, window, 'last')\n",
    "    cows_pred_VAR, calves_pred_VAR = rolling_predictions(endog_diff, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, window, 'VAR')\n",
    "    \n",
    "    pred_df[f'cows_pred_last_{window}'] = cows_pred_last\n",
    "    pred_df[f'calves_pred_last_{window}'] = calves_pred_last\n",
    "    pred_df[f'cows_pred_VAR_{window}'] = cows_pred_VAR\n",
    "    pred_df[f'calves_pred_VAR_{window}'] = calves_pred_VAR\n",
    "    \n",
    "pred_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae500f1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    \n",
    "    ax.plot(df['cows'], color='blue')\n",
    "    ax.plot(df['calves'], color='orange')\n",
    "    ax.plot(pred_df['cows'], color='blue', label='cows')\n",
    "    ax.plot(pred_df['calves'], color='orange', label='calves')\n",
    "    ax.plot(pred_df[f'cows_pred_last_{windows[i]}'], 'r-.', label='last')\n",
    "    ax.plot(pred_df[f'calves_pred_last_{windows[i]}'], 'r-.')\n",
    "    ax.plot(pred_df[f'cows_pred_VAR_{windows[i]}'], 'k--', label='VAR')\n",
    "    ax.plot(pred_df[f'calves_pred_VAR_{windows[i]}'], 'k--')\n",
    "    ax.legend(loc=2)\n",
    "    ax.set_xlabel('Date')\n",
    "    ax.set_ylabel('$/hundredweight')\n",
    "    ax.axvspan(322, 369, color='#808080', alpha=0.2)\n",
    "    ax.set_xlim(300, 369)\n",
    "    ax.set_title(f'Horizon = {windows[i]}')\n",
    "    \n",
    "plt.setp(axes, xticks=np.arange(300, 369, 12), xticklabels=np.arange(2017, 2023, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4cb8f36",
   "metadata": {},
   "source": [
    "### Evaluation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f29bceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "mae_cow_last_1 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_1'])\n",
    "mae_calves_last_1 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_1'])\n",
    "mae_cow_VAR_1 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VAR_1'])\n",
    "mae_calves_VAR_1 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VAR_1'])\n",
    "\n",
    "mae_cow_last_4 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_4'])\n",
    "mae_calves_last_4 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_4'])\n",
    "mae_cow_VAR_4 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VAR_4'])\n",
    "mae_calves_VAR_4 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VAR_4'])\n",
    "\n",
    "mae_cow_last_6 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_6'])\n",
    "mae_calves_last_6 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_6'])\n",
    "mae_cow_VAR_6 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VAR_6'])\n",
    "mae_calves_VAR_6 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VAR_6'])\n",
    "\n",
    "mae_cow_last_12 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_12'])\n",
    "mae_calves_last_12 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_12'])\n",
    "mae_cow_VAR_12 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VAR_12'])\n",
    "mae_calves_VAR_12 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VAR_12'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0714c64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mae_cows_last = [mae_cow_last_1, mae_cow_last_4, mae_cow_last_6, mae_cow_last_12]\n",
    "mae_calves_last = [mae_calves_last_1, mae_calves_last_4, mae_calves_last_6, mae_calves_last_12]\n",
    "mae_cows_VAR = [mae_cow_VAR_1, mae_cow_VAR_4, mae_cow_VAR_6, mae_cow_VAR_12]\n",
    "mae_calves_VAR = [mae_calves_VAR_1, mae_calves_VAR_4, mae_calves_VAR_6, mae_calves_VAR_12]\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "width = 0.3\n",
    "labels = ['Cows', 'Calves']\n",
    "x = np.arange(len(labels))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    y_last = [mae_cows_last[i], mae_calves_last[i]]\n",
    "    y_VAR = [mae_cows_VAR[i], mae_calves_VAR[i]]\n",
    "    \n",
    "    ax.bar(x - width/2, y_last, width, label='Last')\n",
    "    ax.bar(x + width/2, y_VAR, width, label='VAR')\n",
    "    \n",
    "    ax.set_ylim(0, max(y_last) + 10)\n",
    "    ax.set_ylabel('MAE')\n",
    "    ax.set_xlabel('Target')\n",
    "    ax.set_title(f'Horizon = {windows[i]}')\n",
    "    ax.set_xticks(x)\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.legend(loc='best')\n",
    "    \n",
    "    for i, v in enumerate(y_last):\n",
    "        ax.text(x=i - width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "    for i, v in enumerate(y_VAR):\n",
    "        ax.text(x=i + width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efe4c7d7",
   "metadata": {},
   "source": [
    "## VARMA(p,q) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2784f542",
   "metadata": {},
   "outputs": [],
   "source": [
    "def VARMA_gridsearch(endog, min_p, max_p, min_q, max_q):\n",
    "    \n",
    "    all_p = range(min_p, max_p+1, 1)\n",
    "    # Define the range of values for q\n",
    "\n",
    "    # Make a list of all possible (p,q) combinations\n",
    "    \n",
    "    results = []\n",
    "    \n",
    "    print(f'Fitting {len(all_orders)} unique models')\n",
    "    \n",
    "    for order in tqdm_notebook(all_orders):\n",
    "        try:\n",
    "            # Fit a VARMA model\n",
    "        except:\n",
    "            continue\n",
    "    \n",
    "        results.append([order, model.aic])\n",
    "        \n",
    "    result_df = pd.DataFrame(results)\n",
    "    result_df.columns = ['(p,q)', 'AIC']\n",
    "    \n",
    "    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)\n",
    "    \n",
    "    return result_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8b579f7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ad8fa54b",
   "metadata": {},
   "source": [
    "We don't need to redo the Granger causality test, since $p=6$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1c21d0b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "db0595af",
   "metadata": {},
   "source": [
    "### Residuals analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc6b925",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cows\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef1191ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cf5d664",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calves\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4778a0a1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "12d93513",
   "metadata": {},
   "source": [
    "### Forecasting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa126c4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rolling_predictions(df, last_train_value, train_len, horizon, window, method):\n",
    "    \n",
    "    total_len = train_len + horizon\n",
    "    \n",
    "    if method == 'VARMA':\n",
    "\n",
    "        cows_pred_VARMA = []\n",
    "        calves_pred_VARMA = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            # Initialize the VARMA(p,q) model\n",
    "            \n",
    "            res = model.fit(disp=False)\n",
    "            predictions = res.get_prediction(0, i + window - 1)\n",
    "            \n",
    "            oos_pred_cows = predictions.predicted_mean.iloc[-window:]['cows']\n",
    "            oos_pred_calves = predictions.predicted_mean.iloc[-window:]['calves']\n",
    "            \n",
    "            cows_pred_VARMA.extend(oos_pred_cows)\n",
    "            calves_pred_VARMA.extend(oos_pred_calves)\n",
    "            \n",
    "        cows_pred_VARMA = np.insert(cows_pred_VARMA, 0, last_train_value['cows'])\n",
    "        cows_pred_VARMA = cows_pred_VARMA.cumsum()\n",
    "        \n",
    "        calves_pred_VARMA = np.insert(calves_pred_VARMA, 0, last_train_value['calves'])\n",
    "        calves_pred_VARMA = calves_pred_VARMA.cumsum()\n",
    "        \n",
    "        return cows_pred_VARMA[:horizon], calves_pred_VARMA[:horizon]\n",
    "    \n",
    "    elif method == 'last':\n",
    "        cows_pred_last = []\n",
    "        calves_pred_last = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            \n",
    "            cows_last = df[:i].iloc[-1]['cows']\n",
    "            calves_last = df[:i].iloc[-1]['calves']\n",
    "            \n",
    "            cows_pred_last.extend(cows_last for _ in range(window))\n",
    "            calves_pred_last.extend(calves_last for _ in range(window))\n",
    "        \n",
    "        cows_pred_last = np.insert(cows_pred_last, 0, last_train_value['cows'])\n",
    "        cows_pred_last = cows_pred_last.cumsum()\n",
    "        \n",
    "        calves_pred_last = np.insert(calves_pred_last, 0, last_train_value['calves'])\n",
    "        calves_pred_last = calves_pred_last.cumsum()\n",
    "            \n",
    "        return cows_pred_last[:horizon], calves_pred_last[:horizon]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "302e8a7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_df = df[['cows', 'calves']][-48:]\n",
    "\n",
    "TRAIN_LEN = len(train)\n",
    "HORIZON = len(test)\n",
    "LAST_TRAIN_VALUE = endog.iloc[321]\n",
    "\n",
    "windows = [1, 4, 6, 12]\n",
    "\n",
    "for window in windows:\n",
    "    cows_pred_last, calves_pred_last = rolling_predictions(endog_diff, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, window, 'last')\n",
    "    cows_pred_VARMA, calves_pred_VARMA = rolling_predictions(endog_diff, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, window, 'VARMA')\n",
    "    \n",
    "    pred_df[f'cows_pred_last_{window}'] = cows_pred_last\n",
    "    pred_df[f'calves_pred_last_{window}'] = calves_pred_last\n",
    "    pred_df[f'cows_pred_VARMA_{window}'] = cows_pred_VARMA\n",
    "    pred_df[f'calves_pred_VARMA_{window}'] = calves_pred_VARMA\n",
    "    \n",
    "pred_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dbc2cbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    \n",
    "    ax.plot(df['cows'], color='blue')\n",
    "    ax.plot(df['calves'], color='orange')\n",
    "    ax.plot(pred_df['cows'], color='blue', label='cows')\n",
    "    ax.plot(pred_df['calves'], color='orange', label='calves')\n",
    "    ax.plot(pred_df[f'cows_pred_last_{windows[i]}'], 'r-.', label='last')\n",
    "    ax.plot(pred_df[f'calves_pred_last_{windows[i]}'], 'r-.')\n",
    "    ax.plot(pred_df[f'cows_pred_VARMA_{windows[i]}'], 'k--', label='VARMA')\n",
    "    ax.plot(pred_df[f'calves_pred_VARMA_{windows[i]}'], 'k--')\n",
    "    ax.legend(loc=2)\n",
    "    ax.set_xlabel('Date')\n",
    "    ax.set_ylabel('$/hundredweight')\n",
    "    ax.axvspan(322, 369, color='#808080', alpha=0.2)\n",
    "    ax.set_xlim(300, 369)\n",
    "    ax.set_title(f'Horizon = {windows[i]}')\n",
    "    \n",
    "plt.setp(axes, xticks=np.arange(300, 369, 12), xticklabels=np.arange(2017, 2023, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2618ad33",
   "metadata": {},
   "source": [
    "### Evaluation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f59e12d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "mae_cow_last_1 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_1'])\n",
    "mae_calves_last_1 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_1'])\n",
    "mae_cow_VARMA_1 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VARMA_1'])\n",
    "mae_calves_VARMA_1 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VARMA_1'])\n",
    "\n",
    "mae_cow_last_4 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_4'])\n",
    "mae_calves_last_4 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_4'])\n",
    "mae_cow_VARMA_4 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VARMA_4'])\n",
    "mae_calves_VARMA_4 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VARMA_4'])\n",
    "\n",
    "mae_cow_last_6 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_6'])\n",
    "mae_calves_last_6 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_6'])\n",
    "mae_cow_VARMA_6 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VARMA_6'])\n",
    "mae_calves_VARMA_6 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VARMA_6'])\n",
    "\n",
    "mae_cow_last_12 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last_12'])\n",
    "mae_calves_last_12 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last_12'])\n",
    "mae_cow_VARMA_12 = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VARMA_12'])\n",
    "mae_calves_VARMA_12 = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VARMA_12'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec8032ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "mae_cows_last = [mae_cow_last_1, mae_cow_last_4, mae_cow_last_6, mae_cow_last_12]\n",
    "mae_calves_last = [mae_calves_last_1, mae_calves_last_4, mae_calves_last_6, mae_calves_last_12]\n",
    "mae_cows_VARMA = [mae_cow_VARMA_1, mae_cow_VARMA_4, mae_cow_VARMA_6, mae_cow_VARMA_12]\n",
    "mae_calves_VARMA = [mae_calves_VARMA_1, mae_calves_VARMA_4, mae_calves_VARMA_6, mae_calves_VARMA_12]\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "width = 0.3\n",
    "labels = ['Cows', 'Calves']\n",
    "x = np.arange(len(labels))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    y_last = [mae_cows_last[i], mae_calves_last[i]]\n",
    "    y_VARMA = [mae_cows_VARMA[i], mae_calves_VARMA[i]]\n",
    "    \n",
    "    ax.bar(x - width/2, y_last, width, label='Last')\n",
    "    ax.bar(x + width/2, y_VARMA, width, label='VARMA')\n",
    "    \n",
    "    ax.set_ylim(0, max(y_last) + 10)\n",
    "    ax.set_ylabel('MAE')\n",
    "    ax.set_xlabel('Target')\n",
    "    ax.set_title(f'Horizon = {windows[i]}')\n",
    "    ax.set_xticks(x)\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.legend(loc='best')\n",
    "    \n",
    "    for i, v in enumerate(y_last):\n",
    "        ax.text(x=i - width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "    for i, v in enumerate(y_VARMA):\n",
    "        ax.text(x=i + width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75b6c788",
   "metadata": {},
   "source": [
    "## VARMAX "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aa3f100",
   "metadata": {},
   "outputs": [],
   "source": [
    "def VARMAX_gridsearch(endog, exog, min_p, max_p, min_q, max_q):\n",
    "    \n",
    "    all_p = range(min_p, max_p+1, 1)\n",
    "    all_q = range(min_q, max_q+1, 1)\n",
    "    all_orders = list(product(all_p, all_q))\n",
    "    \n",
    "    print(f'Fitting {len(all_orders)} unique models')\n",
    "    \n",
    "    results = []\n",
    "    \n",
    "    for order in tqdm_notebook(all_orders):\n",
    "        try:\n",
    "            # Fit a VARMAX model\n",
    "            \n",
    "        except:\n",
    "            continue\n",
    "    \n",
    "        results.append([order, model.aic])\n",
    "        \n",
    "    result_df = pd.DataFrame(results)\n",
    "    result_df.columns = ['(p,q)', 'AIC']\n",
    "    \n",
    "    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)\n",
    "    \n",
    "    return result_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5288fbd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get your exogenous variable\n",
    "# Index of exog must align with the training set\n",
    "\n",
    "# Run the grid search\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49667c05",
   "metadata": {},
   "source": [
    "We don't need to redo the Granger causality test since $p=6$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e143410",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ade7557e",
   "metadata": {},
   "source": [
    "### Residuals analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "118eea4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cows\n",
    "VARMAX_model_fit.plot_diagnostics(variable=1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b129865",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0178b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calves\n",
    "VARMAX_model_fit.plot_diagnostics(variable=1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13b59a54",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9ef74086",
   "metadata": {},
   "source": [
    "### Forecasting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a848deb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rolling_predictions(df, exog, last_train_value, train_len, horizon, window, method):\n",
    "    \n",
    "    total_len = train_len + horizon\n",
    "    \n",
    "    if method == 'VARMAX':\n",
    "\n",
    "        cows_pred_VARMAX = []\n",
    "        calves_pred_VARMAX = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            # Get he predictions from VARAMX\n",
    "            \n",
    "            oos_pred_cows = predictions.predicted_mean.iloc[-window:]['cows']\n",
    "            oos_pred_calves = predictions.predicted_mean.iloc[-window:]['calves']\n",
    "            \n",
    "            cows_pred_VARMAX.extend(oos_pred_cows)\n",
    "            calves_pred_VARMAX.extend(oos_pred_calves)\n",
    "            \n",
    "        cows_pred_VARMAX = np.insert(cows_pred_VARMAX, 0, last_train_value['cows'])\n",
    "        cows_pred_VARMAX = cows_pred_VARMAX.cumsum()\n",
    "        \n",
    "        calves_pred_VARMAX = np.insert(calves_pred_VARMAX, 0, last_train_value['calves'])\n",
    "        calves_pred_VARMAX = calves_pred_VARMAX.cumsum()\n",
    "        \n",
    "        return cows_pred_VARMAX[:horizon], calves_pred_VARMAX[:horizon]\n",
    "    \n",
    "    elif method == 'last':\n",
    "        cows_pred_last = []\n",
    "        calves_pred_last = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            \n",
    "            cows_last = df[:i].iloc[-1]['cows']\n",
    "            calves_last = df[:i].iloc[-1]['calves']\n",
    "            \n",
    "            cows_pred_last.extend(cows_last for _ in range(window))\n",
    "            calves_pred_last.extend(calves_last for _ in range(window))\n",
    "        \n",
    "        cows_pred_last = np.insert(cows_pred_last, 0, last_train_value['cows'])\n",
    "        cows_pred_last = cows_pred_last.cumsum()\n",
    "        \n",
    "        calves_pred_last = np.insert(calves_pred_last, 0, last_train_value['calves'])\n",
    "        calves_pred_last = calves_pred_last.cumsum()\n",
    "            \n",
    "        return cows_pred_last[:horizon], calves_pred_last[:horizon]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d48c52",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_df = df[['cows', 'calves']][-48:]\n",
    "exog = df['barley'][1:]\n",
    "\n",
    "TRAIN_LEN = len(train)\n",
    "HORIZON = len(test)\n",
    "LAST_TRAIN_VALUE = endog.iloc[321]\n",
    "WINDOW = 1\n",
    "\n",
    "cows_pred_last, calves_pred_last = rolling_predictions(endog_diff, exog, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, WINDOW, 'last')\n",
    "cows_pred_VARMAX, calves_pred_VARMAX = rolling_predictions(endog_diff, exog, LAST_TRAIN_VALUE, TRAIN_LEN, HORIZON, WINDOW, 'VARMAX')\n",
    "    \n",
    "pred_df[f'cows_pred_last'] = cows_pred_last\n",
    "pred_df[f'calves_pred_last'] = calves_pred_last\n",
    "pred_df[f'cows_pred_VARMAX'] = cows_pred_VARMAX\n",
    "pred_df[f'calves_pred_VARMAX'] = calves_pred_VARMAX\n",
    "    \n",
    "pred_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14c92cd2",
   "metadata": {},
   "source": [
    "### Evaluation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebd835f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "mae_cow_last = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_last'])\n",
    "mae_calves_last = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_last'])\n",
    "mae_cow_VARMAX = mean_absolute_error(pred_df['cows'], pred_df['cows_pred_VARMAX'])\n",
    "mae_calves_VARMAX = mean_absolute_error(pred_df['calves'], pred_df['calves_pred_VARMAX'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8e7ef68",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baade172",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "width = 0.3\n",
    "labels = ['Cows', 'Calves']\n",
    "x = np.arange(len(labels))\n",
    "\n",
    "y_last = [mae_cow_last, mae_calves_last]\n",
    "y_VARMAX = [mae_cow_VARMAX, mae_calves_VARMAX]\n",
    "\n",
    "ax.bar(x - width/2, y_last, width, label='Last')\n",
    "ax.bar(x + width/2, y_VARMAX, width, label='VARMAX')\n",
    "\n",
    "ax.set_ylim(0, 15)\n",
    "ax.set_ylabel('MAE')\n",
    "ax.set_xlabel('Target')\n",
    "ax.set_title(f'Horizon = 1')\n",
    "ax.set_xticks(x)\n",
    "ax.set_xticklabels(labels)\n",
    "ax.legend(loc='best')\n",
    "\n",
    "for i, v in enumerate(y_last):\n",
    "    ax.text(x=i - width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "for i, v in enumerate(y_VARMAX):\n",
    "    ax.text(x=i + width/2, y=v+1, s=str(round(v, 2)), ha='center')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f3892b2",
   "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.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}