08_Forecasting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/ascarrambad/ml-21-22/blob/main/08_Forecasting.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7_JbbGbKA-0s"
   },
   "source": [
    "# Machine Learning SP 2021/2022\n",
    "\n",
    "Prof. Cesare Alippi<br>\n",
    "Giorgia Adorni ([`giorgia.adorni@usi.ch`](mailto:giorgia.adorni@usi.ch))<br>\n",
    "Luca Butera ([`luca.butera@usi.ch`](mailto:luca.butera@usi.ch))<br>\n",
    "Matteo Riva ([`matteo.riva@usi.ch`](mailto:matteo.riva@usi.ch))\n",
    "\n",
    "---\n",
    "# Lab 08: Forecasting\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vOKpKrBLq79U"
   },
   "source": [
    "## Forecasting sunspots\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/4/47/Solar_Archipelago_-_Flickr_-_NASA_Goddard_Photo_and_Video.jpg)\n",
    "\n",
    "As a case study we will cosider the problem of forecasting the number of **sunspots** using a dataset collected by the solar physics research department of the Royal Observatory of Belgium ([link](http://sidc.oma.be/)).\n",
    "\n",
    "According to Wikipedia:\n",
    "\n",
    "> Sunspots are temporary phenomena on the Sun's photosphere that appear as spots darker than the surrounding areas. They are regions of reduced surface temperature caused by concentrations of magnetic field flux that inhibit convection. Sunspots usually appear in pairs of opposite magnetic polarity. Their number varies according to the approximately 11-year solar cycle.\n",
    "\n",
    "The datset is avaiable on Kaggle ([link](https://www.kaggle.com/robervalt/sunspots)) and frequently updated. \n",
    "\n",
    "For semplicity (and reproducibility), we loaded a snapshot of the dataset in the repo of the course.\n",
    "\n",
    "We start by installing the last version of statsmodels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ePdaVOZ7HvWq"
   },
   "outputs": [],
   "source": [
    "%pip install statsmodels -U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jqepj7EelDdS"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "url = \"https://raw.githubusercontent.com/ascarrambad/ml-21-22/main/data/sunspots.csv\"\n",
    "\n",
    "df = pd.read_csv(url, \n",
    "                 parse_dates=True, \n",
    "                 index_col='Date', \n",
    "                 usecols=['Date', 'Monthly Mean Total Sunspot Number'])\n",
    "df = df.rename(columns={'Monthly Mean Total Sunspot Number':'sunspots'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "77JPRz5j5v5d"
   },
   "source": [
    "Let's have a look at the data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4HtmxuqnwgAc"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,5))\n",
    "plt.plot(df.index, df['sunspots'])\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('n. sunspots')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8Gpm5wdb51dR"
   },
   "source": [
    "First of all let's split our data for model evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "N_i8YqCu51Oc"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# NB: DO NOT SHUFFLE THE DATA!!\n",
    "train_df, test_df = train_test_split(df, test_size=0.2, shuffle=False)\n",
    "\n",
    "# scale the data\n",
    "\n",
    "mean = train_df.values.mean()\n",
    "std = train_df.values.std()\n",
    "\n",
    "scale = lambda x : (x - mean) / std\n",
    "inv_scale = lambda x : x * std + mean\n",
    "\n",
    "train_df = scale(train_df)\n",
    "test_df = scale(test_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZTc5gubv2kEc"
   },
   "source": [
    "## A simple baseline\n",
    "\n",
    "Simply use the value at the previous time-step to predict the next one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "iWMByMb-2r4h"
   },
   "outputs": [],
   "source": [
    "actuals = inv_scale(test_df.values.ravel())\n",
    "\n",
    "mse = np.square(actuals[:-1] - actuals[1:]).mean()\n",
    "print(f\"BASELINE MSE: {mse}\")\n",
    "mae = np.abs(actuals[:-1] - actuals[1:]).mean()\n",
    "print(f\"BASELINE MAE: {mae}\")\n",
    "\n",
    "\n",
    "plt.figure(figsize=(20, 5))\n",
    "plt.plot(actuals[1:], label='actual')\n",
    "plt.plot(actuals[:-1], label='pred')\n",
    "plt.xlabel('idx')\n",
    "plt.ylabel('n. sunspots')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jnwCSWlr67Us"
   },
   "source": [
    "### Linear models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Fkc_AFxk6-Vn"
   },
   "source": [
    "Let's start with an AR model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "YDBTC18BSHdl"
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima.model import ARIMA\n",
    "\n",
    "#ARIMA(p,i,q)\n",
    "model = ARIMA(train_df, order=(3, 0, 0), trend='n')\n",
    "\n",
    "fit = model.fit()\n",
    "print(fit.summary())\n",
    "\n",
    "# appends test_df to avaiable data\n",
    "fit = fit.append(test_df, refit=False)\n",
    "\n",
    "actuals = test_df.values.ravel()\n",
    "preds = fit.predict(start=len(train_df)).values.ravel()\n",
    "\n",
    "mse = np.square(inv_scale(preds) - inv_scale(actuals)).mean()\n",
    "print(f\"TEST MSE: {mse}\")\n",
    "mae = np.abs(inv_scale(preds) - inv_scale(actuals)).mean()\n",
    "print(f\"TEST MAE: {mae}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ELKxpf-JCZaB"
   },
   "outputs": [],
   "source": [
    "preds = fit.predict(start=len(train_df), dynamic=False)\n",
    "\n",
    "plt.figure(figsize=(20, 5))\n",
    "plt.plot(inv_scale(test_df), label='actual')\n",
    "plt.plot(inv_scale(preds), label='pred')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('n. sunspots')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Y1wk7XFWFIWB"
   },
   "outputs": [],
   "source": [
    "# one-step prediction\n",
    "y_pred = fit.predict(start=len(train_df), dynamic=False)\n",
    "\n",
    "# recursie prediction\n",
    "n_last = 300\n",
    "y_new = fit.predict(start=len(train_df) + len(test_df) - n_last, dynamic=True)\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(test_df, label=\"actual\")\n",
    "plt.plot(y_pred, label=\"pred\")\n",
    "plt.plot(y_new, label=\"new\")\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('n. sunspots')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XzmLzbLnklWv"
   },
   "source": [
    "## Nonlinear Autoregressive models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wFNM5ofZ5QNQ"
   },
   "source": [
    "### Let's divide data into windows\n",
    "\n",
    "```\n",
    "            historical data                   | future observations\n",
    "--------------------------------------------------------------------\n",
    "x(1) x(2) ... x(t-p) x(t-p+1) ... x(t-1) x(t) | x(t+1)  x(t+2) ...\n",
    "                    \\________________________/| \\____/\n",
    "                       time window            | value to \n",
    "                                              | be predicted\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "c7CmFQr1rzEg"
   },
   "outputs": [],
   "source": [
    "# Time windows\n",
    "def get_time_windows(sequence, window_size):\n",
    "    time = np.arange(sequence.shape[0])\n",
    "    xseq = []\n",
    "    yseq = []\n",
    "    for t in time[:-window_size]:\n",
    "        xseq.append(sequence[t:t+window_size])\n",
    "        yseq.append(sequence[t+window_size])\n",
    "    xseq = np.array(xseq)\n",
    "    yseq = np.array(yseq)\n",
    "    #train-test split\n",
    "    return xseq, yseq.ravel()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "MBK7me3t0xWi"
   },
   "outputs": [],
   "source": [
    "# Prepare the data\n",
    "p = 24\n",
    "\n",
    "x_train, y_train = get_time_windows(sequence=train_df.values.ravel(), window_size=p)\n",
    "x_test, y_test = get_time_windows(sequence=test_df.values.ravel(), window_size=p)\n",
    "print(\"training data\", x_train.shape, y_train.shape)\n",
    "print(\"test data\", x_test.shape, y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "y3MzRgXqKW2a"
   },
   "source": [
    "We can use the usual syntax to build a neural network with keras and train it to predict the next value in the time-series.\n",
    "\n",
    "We can formulate it as a regression problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EvwZ_UHLsk5S"
   },
   "outputs": [],
   "source": [
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Reshape\n",
    "\n",
    "# Build a model \n",
    "def get_model(win_size):\n",
    "    model = Sequential()\n",
    "    model.add(Dense(64, activation=\"relu\", input_shape=(win_size,)))\n",
    "    model.add(Dense(1, activation=\"linear\"))\n",
    "\n",
    "    model.compile(optimizer=\"adam\", loss=\"mse\")\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aqq8ZtvLtKL1"
   },
   "outputs": [],
   "source": [
    "# Prepare the data\n",
    "p = 24\n",
    "\n",
    "x_train, y_train = get_time_windows(sequence=train_df.values.ravel(), window_size=p)\n",
    "x_test, y_test = get_time_windows(sequence=test_df.values.ravel(), window_size=p)\n",
    "\n",
    "# Train the model\n",
    "model = get_model(p)\n",
    "model.fit(x_train, y_train, epochs=50)\n",
    "\n",
    "preds = model.predict(x_test).ravel()\n",
    "\n",
    "mse = np.square(inv_scale(preds) - inv_scale(y_test)).mean()\n",
    "print(f\"TEST MSE: {mse}\")\n",
    "mae = np.abs(inv_scale(preds) - inv_scale(y_test)).mean()\n",
    "print(f\"TEST MAE: {mae}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kEFtODpIuXor"
   },
   "outputs": [],
   "source": [
    "# Predictions\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(np.arange(y_test.shape[0]), y_test, label=\"actual\")\n",
    "plt.plot(np.arange(y_test.shape[0]), y_pred, label=\"pred\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n499rYGT61rm"
   },
   "source": [
    "It looks like that the NAR model is not very good.\n",
    "\n",
    "Let's see what happens with recursive predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "IRzxDWsgusvt"
   },
   "outputs": [],
   "source": [
    "# Predict on our prediction\n",
    "def predict_rec(n, model, x_init):\n",
    "    p = x_init.shape[0]\n",
    "\n",
    "    y_new = np.zeros(p+n)\n",
    "    y_new[:p] = x_init.copy()\n",
    "\n",
    "    for j in range(n):\n",
    "        y_tmp = model.predict(y_new[j:j+p][None, ...])\n",
    "        y_new[j+p] = y_tmp[-1,0]\n",
    "    return y_new[p:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DO3CeC-D6zyv"
   },
   "outputs": [],
   "source": [
    "# one-step prediction\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "# one-step prediction\n",
    "n_last = 300\n",
    "y_new = predict_rec(n=n_last, model=model, x_init=x_test[-n_last-1])\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(np.arange(y_test.shape[0]), y_test, label=\"ground truth\")\n",
    "plt.plot(np.arange(y_test.shape[0]), y_pred, label=\"predicted\")\n",
    "plt.plot(np.arange(y_test.shape[0]-n_last, y_test.shape[0]), y_new, label=\"new\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QFRE6_e5HB5P"
   },
   "source": [
    "## Recurrent neural networks\n",
    "\n",
    "Finally let's try out a more advanced architecture: LSTM networks.\n",
    "\n",
    "Implementing LSTMs is again easy with keras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KxsN2Tc1h1YB"
   },
   "outputs": [],
   "source": [
    "from tensorflow.keras.layers import LSTM, Reshape\n",
    "\n",
    "# Build a model \n",
    "def get_model(win_size):\n",
    "    model = Sequential()\n",
    "    model.add(Reshape((-1, 1), input_shape=(None,))) # [batch, time, features]\n",
    "    model.add(LSTM(16))\n",
    "    model.add(Dense(1, activation=\"linear\"))\n",
    "\n",
    "    model.compile(optimizer=\"adam\", loss=\"mse\")\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-YIOYFNqnv1k"
   },
   "outputs": [],
   "source": [
    "# Prepare the data\n",
    "p = 150\n",
    "\n",
    "x_train, y_train = get_time_windows(sequence=train_df.values.ravel(), window_size=p)\n",
    "x_test, y_test = get_time_windows(sequence=test_df.values.ravel(), window_size=p)\n",
    "\n",
    "# Train the model\n",
    "model = get_model(p)\n",
    "model.fit(x_train, y_train, epochs=50)\n",
    "\n",
    "preds = model.predict(x_test).ravel()\n",
    "\n",
    "mse = np.square(inv_scale(preds) - inv_scale(y_test)).mean()\n",
    "print(f\"TEST MSE: {mse}\")\n",
    "\n",
    "mae = np.abs(inv_scale(preds) - inv_scale(y_test)).mean()\n",
    "print(f\"TEST MSE: {mae}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZnmplRAnpLzv"
   },
   "outputs": [],
   "source": [
    "# Predictions\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(np.arange(y_test.shape[0]), y_test, label=\"actual\")\n",
    "plt.plot(np.arange(y_test.shape[0]), y_pred, label=\"pred\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1DAcEwEus9Pi"
   },
   "outputs": [],
   "source": [
    "# Predict on our prediction\n",
    "def predict_rec(n, model, x_init):\n",
    "    p = x_init.shape[0]\n",
    "\n",
    "    y_new = np.zeros(p+n)\n",
    "    y_new[:p] = x_init.copy()\n",
    "\n",
    "    for j in range(n):\n",
    "        y_tmp = model.predict(y_new[:j+p][None, ...])  # the rnn can handle sequences of arbitrary lenght\n",
    "        y_new[j+p] = y_tmp[-1,0]\n",
    "    return y_new[p:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OEM-QtPFtzZt"
   },
   "outputs": [],
   "source": [
    "# one-step prediction\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "# one-step prediction\n",
    "n_last = 150\n",
    "y_new = predict_rec(n=n_last, model=model, x_init=x_test[-n_last-1])\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(np.arange(y_test.shape[0]), y_test, label=\"ground truth\")\n",
    "plt.plot(np.arange(y_test.shape[0]), y_pred, label=\"predicted\")\n",
    "plt.plot(np.arange(y_test.shape[0]-n_last, y_test.shape[0]), y_new, label=\"new\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "s9mqVmE1LBjl"
   },
   "source": [
    "Many more things left to explore:\n",
    "\n",
    "* take seasonality into account;\n",
    "* add exogenous variables (e.g., date/time of the year/...);\n",
    "* multi-step forecasting;\n",
    "* try CNNs;\n",
    "* ...\n",
    "\n",
    "Give a look at [prophet](https://facebook.github.io/prophet/)."
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "08_Forecasting.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}