Task-Time_Series_modeling_Beijing_PM2.5_handout.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling air pollution\n",
    "\n",
    "Let us try to predict Beijing's air pollution, especially [PM2.5](https://en.wikipedia.org/wiki/Particulates) values in advance!\n",
    "\n",
    "\n",
    "Inspiration comes from [here](https://machinelearningmastery.com/multivariate-time-series-forecasting-lstms-keras/).\n",
    "\n",
    "## Dataset\n",
    "\n",
    "[Beijing PM2.5 Data Data Set](https://archive.ics.uci.edu/ml/datasets/Beijing+PM2.5+Data)\n",
    "\n",
    "### Columns of the dataset:\n",
    "\n",
    "**No:** row number\n",
    "\n",
    "**year:** year of data in this row\n",
    "\n",
    "**month:** month of data in this row\n",
    "\n",
    "**day:** day of data in this row\n",
    "\n",
    "**hour:** hour of data in this row\n",
    "\n",
    "**pm2.5:** PM2.5 concentration\n",
    "\n",
    "**DEWP:** Dew Point\n",
    "\n",
    "**TEMP:** Temperature\n",
    "\n",
    "**PRES:** Pressure\n",
    "\n",
    "**cbwd:** Combined wind direction\n",
    "\n",
    "**Iws:** Cumulated wind speed\n",
    "\n",
    "**Is:** Cumulated hours of snow\n",
    "\n",
    "**Ir:** Cumulated hours of rain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:24:58.995243Z",
     "start_time": "2019-05-02T15:24:58.033983Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!wget https://raw.githubusercontent.com/jbrownlee/Datasets/master/pollution.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.287408Z",
     "start_time": "2019-05-02T15:24:59.003580Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!pip install seglearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.659548Z",
     "start_time": "2019-05-02T15:25:00.294723Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import csv\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#import warnings\n",
    "#warnings.filterwarnings(\"ignore\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.665793Z",
     "start_time": "2019-05-02T15:25:00.661169Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sniff_format(location):\n",
    "    with open(location, newline='') as csvfile:\n",
    "        sniffer = csv.Sniffer()\n",
    "        sample = csvfile.read(1024)\n",
    "        dialect = sniffer.sniff(sample)\n",
    "        header = sniffer.has_header(sample)\n",
    "        if header:\n",
    "            header=0\n",
    "        else:\n",
    "            header=None\n",
    "\n",
    "    return {\"dialect\":dialect, \"header\":header}\n",
    "\n",
    "def describe_full(df):\n",
    "    #pd.options.display.float_format = '{:.2f}'.format\n",
    "    dtypes_description=pd.DataFrame(dict(df.dtypes),[\"dtypes\"])\n",
    "    na_description = pd.DataFrame(dict(df.isna().sum()),[\"NA-s\"])\n",
    "    na_percent = ((pd.DataFrame(dict(df.isna().sum()),[\"NA%\"])/len(df))*100).round(decimals=2)\n",
    "    description = df.describe(include='all')\n",
    "    full_description = dtypes_description.append(na_description).append(na_percent).append(description).replace(np.nan, '', regex=True)\n",
    "\n",
    "    mask = full_description.loc[\"freq\",:]==1\n",
    "    full_description.at[[\"top\"],mask.index[mask]]=\"\"\n",
    "    #TODO: scientific notation - could be nicer\n",
    "    \n",
    "    return full_description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.835900Z",
     "start_time": "2019-05-02T15:25:00.667466Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "csv_format = sniff_format(\"pollution.csv\")\n",
    "\n",
    "df = pd.read_csv(\"pollution.csv\",header=csv_format[\"header\"],dialect=csv_format[\"dialect\"])\n",
    "\n",
    "#There is a warning that would be worth investigationg, but for now, let's ignore it\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "describe_full(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.845590Z",
     "start_time": "2019-05-02T15:25:00.837373Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Safety assertions to ensure, that No is as it seems an index column sorted, series\n",
    "pd.testing.assert_series_equal(df.No,df.No.sort_values())\n",
    "np.testing.assert_array_equal((df.No-df.No[0]).values, np.arange(0,df.shape[0],1))\n",
    "#Don't really need these if we decide to go for date based indexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:00.992198Z",
     "start_time": "2019-05-02T15:25:00.847885Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df.drop(\"No\",axis=1, inplace=True)\n",
    "df[\"date\"]= pd.to_datetime(df['year'].astype(str)+'-'+df['month'].astype(str)+\"-\"+df[\"day\"].astype(str)+\"T\"+df[\"hour\"].astype(str).apply(lambda x: x.zfill(2)+\":00\"))\n",
    "df.set_index(df.date, inplace=True)\n",
    "df.drop(\"date\", axis=1, inplace=True)\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encoding day of week\n",
    "\n",
    "We explicitly encode the day of week, since we assume that weekends and workdays behave differently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:01.000339Z",
     "start_time": "2019-05-02T15:25:00.993439Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df[\"dayofweek\"]=df.index.dayofweek+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decision about NaN-s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:01.354703Z",
     "start_time": "2019-05-02T15:25:01.001819Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig = plt.gcf()\n",
    "fig.set_size_inches(10,7)\n",
    "df[\"pm2.5\"].isnull().astype(float).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:01.432487Z",
     "start_time": "2019-05-02T15:25:01.356144Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def isnullsum(df):\n",
    "    return df.isnull().sum()\n",
    "print(\"1-----\")\n",
    "print(\"% NaN datapoints per year:\",(df.groupby(df.index.year)[\"pm2.5\"].apply(isnullsum)/df.groupby(df.index.year)[\"pm2.5\"].count())*100.0)\n",
    "\n",
    "print(\"2--------------------------------\")\n",
    "\n",
    "\n",
    "print(\"% NaN datapoints per month:\",(df.groupby(df.index.month)[\"pm2.5\"].apply(isnullsum)/df.groupby(df.index.month)[\"pm2.5\"].count())*100.0)\n",
    "\n",
    "print(\"3--------------------------------\")\n",
    "\n",
    "print(\"% NaN datapoints per dayofweek:\",(df.groupby(df.index.dayofweek)[\"pm2.5\"].apply(isnullsum)/df.groupby(df.index.dayofweek)[\"pm2.5\"].count())*100.0)\n",
    "\n",
    "print(\"4--------------------------------\")\n",
    "print(\"% NaN datapoints per hour:\",(df.groupby(df.index.hour)[\"pm2.5\"].apply(isnullsum)/df.groupby(df.index.hour)[\"pm2.5\"].count())*100.0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After examining the NaN values in pm2.5, we see no obvious temporal pattern. This is cause for worry, since by simply dropping the rows with NaN values, we can destroy the temporal coherence of the data, hence **data imputation is desirable.**\n",
    "\n",
    "The autocorrelation charts below imply, that it is not unreasonable to take the previous value to fill NaN-s (high autocorrelation with the previous timestep)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:01.471480Z",
     "start_time": "2019-05-02T15:25:01.433895Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "df.fillna(method='ffill', inplace=True)\n",
    "\n",
    "print(df.isnull().sum())\n",
    "\n",
    "df.dropna(inplace=True)\n",
    "\n",
    "print(df.isnull().sum())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examining autocorrelations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:25.702761Z",
     "start_time": "2019-05-02T15:25:01.472908Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from statsmodels.graphics.tsaplots import plot_pacf\n",
    "\n",
    "#columns = [] #use this for speedup\n",
    "columns = [\"pm2.5\",\"DEWP\",\"TEMP\",\"PRES\",\"Iws\",\"Is\",\"Ir\"]\n",
    "\n",
    "for col in columns:\n",
    "    plt.figure()\n",
    "    plot_pacf(df[col].dropna(), lags=200, zero=False)\n",
    "    \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What do we see?\n",
    "\n",
    "Well, the fact, that we don't see.\n",
    "\n",
    "Or more precisely: smog (and weather) is slow to move, it is extremely strongly autocorrelated with itself one-two hours before, so in order to at least be able to see some autocorrelation structure beyond this, we need to filter out the first some hours from our autocorrelation analysis.\n",
    "\n",
    "(By the way, that's why we don't stick to the prediction of the next hour as in the original \"inspiration\" blogpost. Would not be too relevant...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:25.715038Z",
     "start_time": "2019-05-02T15:25:25.705295Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from statsmodels.graphics.tsaplots import _prepare_data_corr_plot, _plot_corr\n",
    "import statsmodels.graphics.utils as utils\n",
    "from statsmodels.tsa.stattools import pacf\n",
    "\n",
    "def plot_pacf_drop(x, ax=None, lags=None, alpha=.05, method='ywadjusted',\n",
    "              use_vlines=True, title='Partial Autocorrelation', zero=True,\n",
    "              vlines_kwargs=None, drop_no=0, **kwargs):\n",
    "    \n",
    "    lags_orig=lags\n",
    "    fig, ax = utils.create_mpl_ax(ax)\n",
    "    vlines_kwargs = {} if vlines_kwargs is None else vlines_kwargs\n",
    "    lags, nlags, irregular = _prepare_data_corr_plot(x, lags, zero)\n",
    "    confint = None\n",
    "    if alpha is None:\n",
    "        acf_x = pacf(x, nlags=nlags, alpha=alpha, method=method)\n",
    "    else:\n",
    "        acf_x, confint = pacf(x, nlags=nlags, alpha=alpha, method=method)\n",
    "\n",
    "    if drop_no:\n",
    "        acf_x = acf_x[drop_no+1:]\n",
    "        confint = confint[drop_no+1:]\n",
    "        lags, nlags, irregular = _prepare_data_corr_plot(x, lags_orig-drop_no, zero)\n",
    "        \n",
    "    _plot_corr(ax, title, acf_x, confint, lags, False, use_vlines,\n",
    "               vlines_kwargs, **kwargs)\n",
    "\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:56.406955Z",
     "start_time": "2019-05-02T15:25:25.718313Z"
    },
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#columns = [] #use this for speedup\n",
    "columns = [\"pm2.5\",\"DEWP\",\"TEMP\",\"PRES\",\"Iws\",\"Is\",\"Ir\"]\n",
    "\n",
    "for col in columns:\n",
    "    plt.figure()\n",
    "    plot_pacf_drop(df[col].dropna(), lags=200, drop_no=3, zero=False)\n",
    "    \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Studying even the filtered charts leaves us in doubt about the possible window for modeling (in case of the classical models), so we will keep 100 as the modeling window (nearly two weeks). This is a parameter that is worth empirically studying later on.\n",
    "\n",
    "It is worth mentioning, that the `pacf` charts would definitely change drastically if we would use some differencing. Since down below we decide not to, we keep it as it is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seasonal decomposition and the question of trends"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:57.851934Z",
     "start_time": "2019-05-02T15:25:56.409704Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "from statsmodels.tsa.tsatools import freq_to_period\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "analysis = seasonal_decompose(df[\"pm2.5\"], period=freq_to_period(df.index.inferred_freq))\n",
    "\n",
    "analysis.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, the default setting (infer periods - hourly) is rather uninformative, so it is maybe worth using some domain knowledge here, and use yearly frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:59.931549Z",
     "start_time": "2019-05-02T15:25:57.854582Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analysis = seasonal_decompose(df[\"pm2.5\"], period=24*365)\n",
    "\n",
    "analysis.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do get the first impression, that there is no overarching simple trend, as well as there are non-trivial seasonal patterns. At a later stage we should investigate differencing regimes, but for now, we leave the data as is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train, valid, test split - before normalization\n",
    "\n",
    "Contamination by the normalization values is a distant possibility, but let's stick to paranoid practices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:25:59.936881Z",
     "start_time": "2019-05-02T15:25:59.933655Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "VALID_AND_TEST_SIZE=0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:00.056548Z",
     "start_time": "2019-05-02T15:25:59.940082Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_else, y_train, y_else = train_test_split(df, df[\"pm2.5\"], test_size=VALID_AND_TEST_SIZE*2, shuffle=False)\n",
    "X_valid, X_test, y_valid, y_test = train_test_split(X_else, y_else, test_size=0.5, shuffle=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have used `temporal_split` from `seglearn`, but that would have cast everything to numpy, so it was more convenient this way for now. Using `seglearn` is encouraged - if we would like to go into classical modeling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-28T08:49:44.040034Z",
     "start_time": "2019-04-28T08:49:43.994879Z"
    }
   },
   "source": [
    "# Data normalization\n",
    "\n",
    "Our default assumption is to use Scikit's minmax scaler for easier learning by neural models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-28T08:49:44.040034Z",
     "start_time": "2019-04-28T08:49:43.994879Z"
    }
   },
   "source": [
    "But there are some exceptions:\n",
    "\n",
    "## How to normalize dates?\n",
    "\n",
    "For the year it is more tricky, it is basically an ordinal.\n",
    "Subtracting the first year is nice, but how to handle the normalization to 0,1?\n",
    "\n",
    "We could use 2018 as a max, but **WE WOULD HAVE TO WRITE A BIG CAVEAT MESSAGE FOR DEPLOY PEOPLE!**\n",
    "\n",
    "So it should be something like  `(df.year - (df.year.min())-1)/((df.year.max()-df.year.min())*2)` (-1 is for avoiding zero, making the life of the network more easy...)\n",
    "\n",
    "For now we stick to the minmax scaler (living risky... :-)\n",
    "\n",
    "For month, day, hour default assumption is, scikit's minmax scaler could work, but we will choose a more elaborate solution from [here](https://ianlondon.github.io/blog/encoding-cyclical-features-24hour-time/). This capitalizes on the circular nature of these quasi ordinals.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:00.065112Z",
     "start_time": "2019-05-02T15:26:00.059198Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "import warnings\n",
    "from sklearn.exceptions import DataConversionWarning\n",
    "warnings.filterwarnings(action='ignore', category=DataConversionWarning)\n",
    "# I literally hate when a standard Scikit function throws big bunches of warnings \n",
    "# - though suppressing them is a dangerous practice. Hence this comment. \n",
    "\n",
    "def minmax_scale(df_x,series_y, normalizers=None):\n",
    "    features_to_minmax = [\"year\",\"pm2.5\",\"DEWP\",\"TEMP\",\"PRES\",\"Iws\",\"Is\",\"Ir\"]\n",
    "\n",
    "    if not normalizers:\n",
    "        normalizers = {}\n",
    "\n",
    "    for feat in features_to_minmax:\n",
    "        if feat not in normalizers:\n",
    "            normalizers[feat] = MinMaxScaler()\n",
    "            normalizers[feat].fit(df_x[feat].values.reshape(-1, 1))\n",
    "        \n",
    "        df_x[feat] = normalizers[feat].transform(df_x[feat].values.reshape(-1, 1))\n",
    "\n",
    "    series_y=normalizers[\"pm2.5\"].transform(series_y.values.reshape(-1, 1))\n",
    "\n",
    "    return df_x, series_y, normalizers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:00.964633Z",
     "start_time": "2019-05-02T15:26:00.067037Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "X_train_norm, y_train_norm, normalizers = minmax_scale(X_train, y_train)\n",
    "X_valid_norm, y_valid_norm, _ = minmax_scale(X_valid, y_valid, normalizers=normalizers)\n",
    "X_test_norm, y_test_norm, _ = minmax_scale(X_test, y_test, normalizers=normalizers)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-28T08:49:44.040034Z",
     "start_time": "2019-04-28T08:49:43.994879Z"
    }
   },
   "source": [
    "## Encoding of ordinals\n",
    "\n",
    "The encoding of `cbwd` is interesting, since it is an ordinal again, or better to say not even that, it has a nice circular topology, so we will use the same sin-cos solution.\n",
    "\n",
    "Problem is, that there is a valid \"zero\" value, marked \"cv\" in there. We are tempted to replace that with 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:00.974581Z",
     "start_time": "2019-05-02T15:26:00.966574Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def encode_cyclicals(df_x):\n",
    "    #\"month\",\"day\",\"hour\", \"cdbw\", \"dayofweek\"\n",
    "    \n",
    "    DIRECTIONS = {\"N\":1.0,\"NE\":2.0, \"E\":3.0, \"SE\":4.0, \"S\":5.0, \"SW\":6.0, \"W\":7.0, \"NW\":8.0, \"cv\":np.nan}\n",
    "\n",
    "    df_x['month_sin'] = np.sin(2*np.pi*df_x.month/12)\n",
    "    df_x['month_cos'] = np.cos(2*np.pi*df_x.month/12)\n",
    "    df_x.drop('month', axis=1, inplace=True)\n",
    "    \n",
    "    df_x['day_sin'] = np.sin(2*np.pi*df_x.day/31)\n",
    "    df_x['day_cos'] = np.cos(2*np.pi*df_x.day/31)\n",
    "    df_x.drop('day', axis=1, inplace=True)\n",
    "\n",
    "    df_x['dayofweek_sin'] = np.sin(2*np.pi*df_x.dayofweek/7)\n",
    "    df_x['dayofweek_cos'] = np.cos(2*np.pi*df_x.dayofweek/7)\n",
    "    df_x.drop('dayofweek', axis=1, inplace=True)\n",
    "    \n",
    "    df_x['hour_sin'] = np.sin(2*np.pi*df_x.hour/24)\n",
    "    df_x['hour_cos'] = np.cos(2*np.pi*df_x.hour/24)\n",
    "    df_x.drop('hour', axis=1, inplace=True)\n",
    "    \n",
    "    df_x.replace({'cbwd': DIRECTIONS}, inplace=True)\n",
    "    df_x['cbwd'] = df_x['cbwd'].astype(np.float64) \n",
    "\n",
    "    df_x['cbwd_sin'] = np.sin(2.0*np.pi*df_x.cbwd/8.0)\n",
    "    df_x['cbwd_sin'].replace(np.nan, 0.0, inplace=True) #Let's handle the case with no wind specially\n",
    "    df_x['cbwd_cos'] = np.cos(2.0*np.pi*df_x.cbwd/8.0)\n",
    "    df_x['cbwd_cos'].replace(np.nan, 0.0, inplace=True) #Let's handle the case with no wind specially\n",
    "    df_x.drop('cbwd', axis=1, inplace=True)\n",
    "    \n",
    "    return df_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:03.103910Z",
     "start_time": "2019-05-02T15:26:00.977468Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train_norm = encode_cyclicals(X_train_norm)\n",
    "X_valid_norm = encode_cyclicals(X_valid_norm)\n",
    "X_test_norm = encode_cyclicals(X_test_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:03.157729Z",
     "start_time": "2019-05-02T15:26:03.105840Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:03.164598Z",
     "start_time": "2019-05-02T15:26:03.159269Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Just in case to ensure we did everything right\n",
    "assert all(x==np.float64 for x in list(X_train_norm.dtypes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would be worth checking with some assertions that the manual normalizers work well. Let's leave it to later work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also worth noting, that the normalizers should be saved and used in production."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating target (y) and \"windows\" (X) for modeling\n",
    "\n",
    "By default we use the next 24 hour value of \"pm2.5\" for prediction, that is, I would like to predict what the pm2.5 will be like **at this hour 24 hours from now.**\n",
    "\n",
    "We use the quite handy **seglearn** package for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of computational reasons, we **use the window of 100 hours** to predict. Classical models would have hard time to accommodate substantially (like 5-10x) context windows, LSTM-s would suffer from the challenge of long term memory. After a basic run of modeling the next big challenge would be to investigate PACF structure more and use eg. stateful LSTM modeling to try to accommodate the large \"lookback\".   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:03.169133Z",
     "start_time": "2019-05-02T15:26:03.165982Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TIME_WINDOW=100\n",
    "FORECAST_DISTANCE=24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.135398Z",
     "start_time": "2019-05-02T15:26:03.170483Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from seglearn.transform import FeatureRep, SegmentXYForecast, last\n",
    "\n",
    "segmenter = SegmentXYForecast(width=TIME_WINDOW, step=1, y_func=last, forecast=FORECAST_DISTANCE)\n",
    "\n",
    "X_train_rolled, y_train_rolled,_=segmenter.fit_transform([X_train_norm.values],[y_train_norm.flatten()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.141218Z",
     "start_time": "2019-05-02T15:26:05.137116Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train_rolled[:1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# For non-sequence models \n",
    "\n",
    "We have to \"flatten\" the data to be able to use classical, non-sequence regression models from Scikit.\n",
    "\n",
    "**We only need to do this for X, any transformation of y is unnecessary.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.154555Z",
     "start_time": "2019-05-02T15:26:05.142809Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train_rolled.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.161443Z",
     "start_time": "2019-05-02T15:26:05.156594Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shape = X_train_rolled.shape\n",
    "X_train_flattened = X_train_rolled.reshape(shape[0],shape[1]*shape[2])\n",
    "X_train_flattened.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.405298Z",
     "start_time": "2019-05-02T15:26:05.163304Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_valid_rolled, y_valid_rolled,_=segmenter.fit_transform([X_valid_norm.values],[y_valid_norm.flatten()])\n",
    "\n",
    "shape = X_valid_rolled.shape\n",
    "X_valid_flattened = X_valid_rolled.reshape(shape[0],shape[1]*shape[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation helper\n",
    "\n",
    "Use this function to evaluate your models **on validation data.**\n",
    "\n",
    "This assumes that your model has the `predict()` function, which is true for **Scikit-learn, XGBoost** and **Keras**, so you can can hand over any of those. \n",
    "\n",
    "A special issue by models optimized by iterative methods is to **get the final model**. **Early stopping and / or model save and reload** can help there.  \n",
    "\n",
    "**WARNING: This is just a basic evaluation scheme, more thorough investigation needed in the future!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.414032Z",
     "start_time": "2019-05-02T15:26:05.407874Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from math import sqrt\n",
    "\n",
    "def evaluate_model(model, X_valid, y_valid_true):\n",
    "    predictions = model.predict(X_valid)\n",
    "    rms = sqrt(mean_squared_error(y_valid_true, predictions))\n",
    "    print(\"Root mean squared error on valid:\",rms)\n",
    "    normalized_rms = normalizers[\"pm2.5\"].inverse_transform(np.array([rms]).reshape(1, -1))[0][0]\n",
    "    print(\"Root mean squared error on valid inverse transformed from normalization:\",normalized_rms)\n",
    "    return normalized_rms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical modeling\n",
    "\n",
    "In \"classical\" modeling we assume a multiple regression case, so we **DO NOT USE time series as such, but the \"flat\" versions of the data** as input. Output is the same. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Baseline - DummyPredictor\n",
    "\n",
    "**TASK Create a dummy predictor as a baseline. Use Scikit-learn's builtin capability to do dummy models in regression case. Use the default setting, that is the prediction of the mean value.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.425330Z",
     "start_time": "2019-05-02T15:26:05.416635Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn...\n",
    "\n",
    "dummy_model = ...\n",
    "\n",
    "dummy_model.fit(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.433386Z",
     "start_time": "2019-05-02T15:26:05.428185Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = evaluate_model(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a RandomForest on raw data\n",
    "\n",
    "**TASK: Fit a RandomForest from Scikit. Please be aware, that the number of trees in the model is having a strong influence on training time.** \n",
    "\n",
    "**Suggestion:** use couple of tens of trees, definitely << 100 to be able to wait it out...\n",
    "\n",
    "**Pro tip:** To utilize all the CPU cores, use the right setting of n_jobs. That speeds things up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:26:05.455678Z",
     "start_time": "2019-05-02T15:26:05.436188Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn...\n",
    "\n",
    "N_ESTIMATORS = ???\n",
    "RANDOM_STATE = 452543634"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:31:10.504738Z",
     "start_time": "2019-05-02T15:26:05.457739Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RF_base_model = ...\n",
    "\n",
    "RF_base_model.fit(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:31:10.647413Z",
     "start_time": "2019-05-02T15:31:10.507882Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = evaluate_model(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a RandomForest on feature transformed data\n",
    "\n",
    "**TASK: Since we use `seglearn`, we can try to capitalize on it's functionality to calculate features from the time time series. Use `FeatureRep` from `seglearn` to transform features, fit a RandomForest and hope for the best!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:32:02.592160Z",
     "start_time": "2019-05-02T15:31:10.654979Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RF_feature_model = ...\n",
    "\n",
    "feature_converter = ...\n",
    "\n",
    "RF_feature_model.fit(feature_converter...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation\n",
    "\n",
    "WARNING: This is just a basic evaluation scheme, more thorough investigation needed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:32:04.242528Z",
     "start_time": "2019-05-02T15:32:02.598939Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = evaluate_model(...,feature_converter..., ...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## XGBoost for speedup\n",
    "\n",
    "Use the XGBoost library to fit gradient boosted trees to the problem. They are usually way quicker to learn and many times at least on par with RandomForests, or better. Let's see!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:32:05.234439Z",
     "start_time": "2019-05-02T15:32:04.248920Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import xgboost as xgb\n",
    "# If in trouble, use !pip install xgboost\n",
    "\n",
    "# XGBoost needs it's custom data format to run quickly\n",
    "dmatrix_train = xgb.DMatrix(data=X_train_flattened,label=y_train_rolled)\n",
    "dmatrix_valid = xgb.DMatrix(data=X_valid_flattened,label=y_valid_rolled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:32:13.647736Z",
     "start_time": "2019-05-02T15:32:05.237071Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = {'objective': 'reg:linear', 'eval_metric': 'rmse', 'n_estimators': ??}\n",
    "\n",
    "evallist = [(dmatrix_valid, 'eval'), (dmatrix_train, 'train')]\n",
    "\n",
    "num_round = ?? #Can easily overfit, experiment with it!\n",
    "\n",
    "xg_reg = xgb.train(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:32:13.654770Z",
     "start_time": "2019-05-02T15:32:13.649668Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = evaluate_model(...,dmatrix...,...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building an LSTM model\n",
    "\n",
    "## Modeling assumptions\n",
    "\n",
    "**TASK:** We believe, that the time dependent structure of this dataset is complex, so we try to use LSTM models from Keras. We are not explicitly utilizing **statefulness**, that is a **major area to be investigated later on**.\n",
    "\n",
    "More information on statefulness can be found [here](https://philipperemy.github.io/keras-stateful-lstm/).\n",
    "\n",
    "\n",
    "Fit an LSTM model **on the time series - non-flat - data!**.\n",
    "\n",
    "Use:\n",
    "1. At least 1 LSTM layer\n",
    "2. A dense layer for output - think about activation! This is a regression case!\n",
    "\n",
    "**Very advisable** - but optional - to use Dropout. You can not use it everywhere, though... Experiment!\n",
    "\n",
    "You are allowed to use functional API, but for this **Sequential API is sufficient.**\n",
    "\n",
    "You **can use LeraningRateScheduler** if you like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:36:43.004415Z",
     "start_time": "2019-05-02T15:36:43.001752Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LSTM_CELL_SIZE=??\n",
    "BATCH_SIZE = ??\n",
    "EPOCHS = ??\n",
    "DROPOUT_RATE=??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:43:36.343431Z",
     "start_time": "2019-05-02T15:36:43.007323Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from tensorflow.keras import Sequential\n",
    "from tensorflow.keras.layers import Dense, Dropout, LSTM\n",
    "from tensorflow.keras import backend as be\n",
    "from tensorflow.keras.callbacks import ModelCheckpoint, LearningRateScheduler\n",
    "\n",
    "column_count=len(X_train_???.columns) #Remember,column count before rolling...\n",
    "\n",
    "be.clear_session()\n",
    "\n",
    "# You might very well be needing it!\n",
    "# Remeber to save only what is worth it from validation perspective...\n",
    "# model_saver = ModelCheckpoint(...)\n",
    "\n",
    "# If you need it...\n",
    "#def schedule(epoch, lr):\n",
    "#    ...\n",
    "#    return lr\n",
    "\n",
    "#lr_scheduler = LearningRateScheduler(schedule)\n",
    "\n",
    "# Build your whole LSTM model here!\n",
    "model = ...\n",
    "\n",
    "#For shape remeber, we have a variable defining the \"window\" and the features in the window...\n",
    "\n",
    "model.compile(loss='mean_squared_error', optimizer=???)\n",
    "# Fit on the train data\n",
    "# USE the batch size parameter!\n",
    "# Use validation data - warning, a tuple of stuff!\n",
    "# Epochs as deemed necessary...\n",
    "# You should avoid shuffling the data maybe.\n",
    "# You can use the callbacks for LR schedule or model saving as seems fit.\n",
    "history = model.fit(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:43:36.500553Z",
     "start_time": "2019-05-02T15:43:36.347851Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(history.history['loss'], label='train')\n",
    "plt.plot(history.history['val_loss'], label='test')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:43:44.053319Z",
     "start_time": "2019-05-02T15:43:36.503640Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# You can use the early stopped model OR load it. \n",
    "# For that you have to import the load function...\n",
    "# IF AND ONLY IF loading, it is good practice to throw out the trash from the graph...\n",
    "# be.clear_session()\n",
    "\n",
    "\n",
    "result = evaluate_model(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-02T15:43:44.057870Z",
     "start_time": "2019-05-02T15:43:44.055366Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assert result < 86.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Things that should be improved\n",
    "\n",
    "- More conclusive investigation of PACF for better time window estimate\n",
    "    - It can well be, that long windows do not add that much to the performance\n",
    "- More interesting features for XGBoost (like from [tsfresh](https://tsfresh.readthedocs.io/en/latest/)), since present features are a disaster\n",
    "- MOST IMPORTANT: **More thorough error / prediction analysis!!!**\n",
    "- LSTM with **Custom iterator with stateful model**\n",
    "- Investigation of different loss function (eg. MAE) for training. (And with it, think abut the importance of extreme values: do we think they are outliers? Are they interesing to predict?)\n",
    "- Investigation of \"teacher forcing\" for LSTM-s in Keras (if it makes sense)\n",
    "\n",
    "# Conclusion\n",
    "\n",
    "Even with decent amount of struggle, the \"dummy\" of always using the mean is very appealing, so it seems, this is not that easy of a task 24 hours in advance. Further investigation of classical as well as neural models remains open!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final test\n",
    "\n",
    "We did not use the final test, since our investigations are not concluded yet. Remember: using it once before project \"go live\" is a good practice!"
   ]
  }
 ],
 "metadata": {
  "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.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}