3_multi_step_FF_vector_output.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi step model (vector output approach)\n",
    "\n",
    "Download zipfile from https://www.dropbox.com/s/pqenrr2mcvl0hk9/GEFCom2014.zip?dl=0 and store in the data folder.\n",
    "\n",
    "In this notebook, we demonstrate how to:\n",
    "- prepare time series data for training a RNN forecasting model\n",
    "- get data in the required shape for the keras API\n",
    "- implement a RNN model in keras to predict the next 24 steps ahead (time *t+1* to *t+24*) in the time series. This model uses recent values of temperature and load as the model input. The model will be trained to output a vector, the elements of which are ordered predictions for future time steps.\n",
    "- enable early stopping to reduce the likelihood of model overfitting\n",
    "- evaluate the model on a test dataset\n",
    "\n",
    "The data in this example is taken from the GEFCom2014 forecasting competition<sup>1</sup>. It consists of 3 years of hourly electricity load and temperature values between 2012 and 2014. The task is to forecast future values of electricity load.\n",
    "\n",
    "<sup>1</sup>Tao Hong, Pierre Pinson, Shu Fan, Hamidreza Zareipour, Alberto Troccoli and Rob J. Hyndman, \"Probabilistic energy forecasting: Global Energy Forecasting Competition 2014 and beyond\", International Journal of Forecasting, vol.32, no.3, pp 896-913, July-September, 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "from glob import glob\n",
    "from collections import UserDict\n",
    "import itertools\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "pd.options.display.float_format = '{:,.2f}'.format\n",
    "np.set_printoptions(precision=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run -i common/load_data.py\n",
    "%run -i common/mape.py\n",
    "%run -i common/TimeSeriesTensor.py\n",
    "%run -i common/create_evaluation_df.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data into Pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "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>load</th>\n",
       "      <th>temp</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2012-01-01 00:00:00</th>\n",
       "      <td>2,698.00</td>\n",
       "      <td>32.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-01-01 01:00:00</th>\n",
       "      <td>2,558.00</td>\n",
       "      <td>32.67</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-01-01 02:00:00</th>\n",
       "      <td>2,444.00</td>\n",
       "      <td>30.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-01-01 03:00:00</th>\n",
       "      <td>2,402.00</td>\n",
       "      <td>31.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-01-01 04:00:00</th>\n",
       "      <td>2,403.00</td>\n",
       "      <td>32.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        load  temp\n",
       "2012-01-01 00:00:00 2,698.00 32.00\n",
       "2012-01-01 01:00:00 2,558.00 32.67\n",
       "2012-01-01 02:00:00 2,444.00 30.00\n",
       "2012-01-01 03:00:00 2,402.00 31.00\n",
       "2012-01-01 04:00:00 2,403.00 32.00"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "if not os.path.exists(os.path.join('data', 'energy.csv')):\n",
    "    %run common/extract_data.py\n",
    "energy = load_data()\n",
    "energy.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "valid_start_dt = '2014-09-01 00:00:00'\n",
    "test_start_dt = '2014-11-01 00:00:00'\n",
    "\n",
    "T = 6\n",
    "HORIZON = 24\n",
    "N_EXPERIMENTS = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train = energy.copy()[energy.index < valid_start_dt][['load', 'temp']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "y_scaler = MinMaxScaler()\n",
    "y_scaler.fit(train[['load']])\n",
    "\n",
    "X_scaler = MinMaxScaler()\n",
    "train[['load', 'temp']] = X_scaler.fit_transform(train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implement the RNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\h5py\\__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "from keras.models import Model, Sequential\n",
    "from keras.layers import GRU, Dense\n",
    "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
    "from keras.optimizers import RMSprop\n",
    "from keras import losses, regularizers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fixed parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "EPOCHS = 100\n",
    "ALPHA = 0.01 # regularization coefficient\n",
    "N_EXPERIMENTS = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tunable parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "LATENT_DIMS = [5, 10, 15]\n",
    "BATCH_SIZES = [8, 16, 32]\n",
    "LEARNING_RATES = [0.01, 0.001, 0.0001]\n",
    "ALL_PARAMS = [LATENT_DIMS, BATCH_SIZES, LEARNING_RATES]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T_VALUES = [3, 7, 14]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = [list(enumerate(x)) for x in ALL_PARAMS]\n",
    "grid = list(itertools.product(*parameters))\n",
    "lengths = [len(T_VALUES)]\n",
    "lengths.extend([len(x) for x in ALL_PARAMS])\n",
    "mapes = np.empty(tuple(lengths))\n",
    "st_errs = np.empty_like(mapes) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_input(T):\n",
    "    \n",
    "    tensor_structure = {'X':(range(-T+1, 1), ['load', 'temp'])}\n",
    "    train_inputs = TimeSeriesTensor(train, 'load', HORIZON, tensor_structure)\n",
    "    X_train = train_inputs.dataframe.as_matrix()[:,HORIZON:]\n",
    "    Y_train = train_inputs['target']\n",
    "    \n",
    "    # Construct validation set (keeping T hours from the training set in order to construct initial features)\n",
    "    look_back_dt = dt.datetime.strptime(valid_start_dt, '%Y-%m-%d %H:%M:%S') - dt.timedelta(hours=T-1)\n",
    "    valid = energy.copy()[(energy.index >=look_back_dt) & (energy.index < test_start_dt)][['load', 'temp']]\n",
    "    valid[['load', 'temp']] = X_scaler.transform(valid)\n",
    "    valid_inputs = TimeSeriesTensor(valid, 'load', HORIZON, tensor_structure)\n",
    "    X_valid = valid_inputs.dataframe.as_matrix()[:,HORIZON:]\n",
    "    Y_valid = valid_inputs['target']\n",
    "    \n",
    "    # Construct test set (keeping T hours from the validation set in order to construct initial features)\n",
    "    look_back_dt = dt.datetime.strptime(test_start_dt, '%Y-%m-%d %H:%M:%S') - dt.timedelta(hours=T-1)\n",
    "    test = energy.copy()[test_start_dt:][['load', 'temp']]\n",
    "    test[['load', 'temp']] = X_scaler.transform(test)\n",
    "    test_inputs = TimeSeriesTensor(test, 'load', HORIZON, tensor_structure)\n",
    "    X_test = test_inputs.dataframe.as_matrix()[:,HORIZON:]\n",
    "    \n",
    "    return X_train, Y_train, X_valid, Y_valid, X_test, test_inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_model(LATENT_DIM, LEARNING_RATE, T, ALPHA):\n",
    "    model = Sequential()\n",
    "    model.add(Dense(LATENT_DIM, activation=\"relu\", input_shape=(2*T,), \\\n",
    "                    kernel_regularizer=regularizers.l2(ALPHA), bias_regularizer=regularizers.l2(ALPHA)))\n",
    "    model.add(Dense(HORIZON, kernel_regularizer=regularizers.l2(ALPHA), bias_regularizer=regularizers.l2(ALPHA)))\n",
    "    optimizer = RMSprop(lr=LEARNING_RATE)\n",
    "    model.compile(optimizer=optimizer, loss=losses.mean_absolute_percentage_error)\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1705\n",
      "0.0886\n",
      "Mean MAPE = 0.1295 +/- 0.0290\n",
      "0.0886\n",
      "0.0982\n",
      "Mean MAPE = 0.0934 +/- 0.0034\n",
      "0.0882\n",
      "0.0897\n",
      "Mean MAPE = 0.0889 +/- 0.0005\n",
      "0.1028\n",
      "0.1072\n",
      "Mean MAPE = 0.1050 +/- 0.0016\n",
      "0.0802\n",
      "0.0970\n",
      "Mean MAPE = 0.0886 +/- 0.0059\n",
      "0.0886\n",
      "0.0882\n",
      "Mean MAPE = 0.0884 +/- 0.0002\n",
      "0.1345\n",
      "0.1366\n",
      "Mean MAPE = 0.1356 +/- 0.0008\n",
      "0.0882\n",
      "0.0886\n",
      "Mean MAPE = 0.0884 +/- 0.0001\n",
      "0.0870\n",
      "0.1008\n",
      "Mean MAPE = 0.0939 +/- 0.0049\n",
      "0.0895\n",
      "0.0935\n",
      "Mean MAPE = 0.0915 +/- 0.0014\n",
      "0.0840\n",
      "0.0826\n",
      "Mean MAPE = 0.0833 +/- 0.0005\n",
      "0.0848\n",
      "0.0849\n",
      "Mean MAPE = 0.0848 +/- 0.0000\n",
      "0.1216\n",
      "0.1255\n",
      "Mean MAPE = 0.1236 +/- 0.0014\n",
      "0.0797\n",
      "0.0914\n",
      "Mean MAPE = 0.0856 +/- 0.0041\n",
      "0.0853\n",
      "0.0888\n",
      "Mean MAPE = 0.0871 +/- 0.0013\n",
      "0.1121\n",
      "0.1794\n",
      "Mean MAPE = 0.1458 +/- 0.0238\n",
      "0.0865\n",
      "0.1018\n",
      "Mean MAPE = 0.0941 +/- 0.0054\n",
      "0.0999\n",
      "0.0886\n",
      "Mean MAPE = 0.0943 +/- 0.0040\n",
      "0.0902\n",
      "0.0790\n",
      "Mean MAPE = 0.0846 +/- 0.0040\n",
      "0.0794\n",
      "0.0907\n",
      "Mean MAPE = 0.0851 +/- 0.0040\n",
      "0.0916\n",
      "0.0833\n",
      "Mean MAPE = 0.0875 +/- 0.0029\n",
      "0.1217\n",
      "0.1181\n",
      "Mean MAPE = 0.1199 +/- 0.0013\n",
      "0.0832\n",
      "0.0938\n",
      "Mean MAPE = 0.0885 +/- 0.0037\n",
      "0.0939\n",
      "0.0846\n",
      "Mean MAPE = 0.0892 +/- 0.0033\n",
      "0.1242\n",
      "0.1199\n",
      "Mean MAPE = 0.1221 +/- 0.0015\n",
      "0.0963\n",
      "0.0972\n",
      "Mean MAPE = 0.0968 +/- 0.0003\n",
      "0.0917\n",
      "0.0898\n",
      "Mean MAPE = 0.0908 +/- 0.0006\n"
     ]
    }
   ],
   "source": [
    "for t, T_val in enumerate(T_VALUES):\n",
    "    \n",
    "    X_train, Y_train, X_valid, Y_valid, X_test, test_inputs = create_input(T_val)\n",
    "      \n",
    "    for (i,LATENT_DIM), (j,BATCH_SIZE), (k,LEARNING_RATE) in grid:\n",
    "    \n",
    "        mapes_param = np.empty(N_EXPERIMENTS)\n",
    "        for ii in range(N_EXPERIMENTS):\n",
    "    \n",
    "            # Initialize the model\n",
    "            model = get_model(BATCH_SIZE, LEARNING_RATE, T_val, ALPHA)\n",
    "            earlystop = EarlyStopping(monitor='val_loss', min_delta=0, patience=5)\n",
    "            best_val = ModelCheckpoint('model_{epoch:02d}.h5', save_best_only=True, mode='min', period=1)\n",
    "    \n",
    "            # Train the model\n",
    "            history = model.fit(X_train, Y_train,\n",
    "                                batch_size=BATCH_SIZE,\n",
    "                                epochs=EPOCHS,\n",
    "                                validation_data=(X_valid, Y_valid),\n",
    "                                callbacks=[earlystop, best_val],\n",
    "                                verbose=0)\n",
    "    \n",
    "            # load the model with the smallest MAPE\n",
    "            best_epoch = np.argmin(np.array(history.history['val_loss']))+1\n",
    "            model.load_weights(\"model_{:02d}.h5\".format(best_epoch))\n",
    "    \n",
    "            predictions = model.predict(X_test)\n",
    "    \n",
    "            # Compute MAPE for each forecast horizon\n",
    "            eval_df = create_evaluation_df(predictions, test_inputs, HORIZON, y_scaler)\n",
    "            eval_df['APE'] = (eval_df['prediction'] - eval_df['actual']).abs() / eval_df['actual']\n",
    "    \n",
    "            # Compute MAPE across all predictions\n",
    "            mapes_param[ii] = mape(eval_df['prediction'], eval_df['actual'])\n",
    "            print('{0:.4f}'.format(mapes_param[ii]))\n",
    "    \n",
    "            for f in glob('model_*.h5'):\n",
    "                os.remove(f)\n",
    "            \n",
    "        mapes[t, i,j,k] = np.mean(mapes_param)\n",
    "        st_errs[t, i,j,k] = np.std(mapes_param)/np.sqrt(N_EXPERIMENTS)\n",
    "    \n",
    "        result = 'Mean MAPE = {0:.4f} +/- {1:.4f}'.format(np.mean(mapes_param), np.std(mapes_param)/np.sqrt(N_EXPERIMENTS))\n",
    "        print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 3, 3)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mapes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[[1.30e-001, 9.34e-002, 8.89e-002],\n",
       "         [1.05e-001, 8.86e-002, 8.84e-002],\n",
       "         [1.36e-001, 8.84e-002, 9.39e-002]],\n",
       "\n",
       "        [[9.15e-002, 8.33e-002, 8.48e-002],\n",
       "         [1.24e-001, 8.56e-002, 8.71e-002],\n",
       "         [1.46e-001, 9.41e-002, 9.43e-002]],\n",
       "\n",
       "        [[8.46e-002, 8.51e-002, 8.75e-002],\n",
       "         [1.20e-001, 8.85e-002, 8.92e-002],\n",
       "         [1.22e-001, 9.68e-002, 9.08e-002]]],\n",
       "\n",
       "\n",
       "       [[[6.22e+175, 1.14e+243, 1.34e+179],\n",
       "         [2.43e-154, 3.68e+180, 1.34e+161],\n",
       "         [4.18e+199, 1.22e-152, 9.79e+199]],\n",
       "\n",
       "        [[4.82e+228, 6.01e-154, 1.11e+200],\n",
       "         [5.29e+199, 3.17e+180, 1.94e+227],\n",
       "         [7.09e+194, 5.98e-154, 2.46e+198]],\n",
       "\n",
       "        [[4.16e+156, 4.24e+175, 4.98e+151],\n",
       "         [6.18e+223, 5.98e-154, 9.04e-315],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311]]],\n",
       "\n",
       "\n",
       "       [[[1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311]],\n",
       "\n",
       "        [[1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311]],\n",
       "\n",
       "        [[1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311],\n",
       "         [1.09e-311, 1.09e-311, 1.09e-311]]]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}