02_Data_Science-ARIMA-Baseline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Citibike ML\n",
    "In this example notebook we show a baseline model using ARIMA.  This provides a starting point for comparison."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Load the Credentials\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dags.snowpark_connection import snowpark_connect\n",
    "session, state_dict = snowpark_connect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import snowflake.snowpark as snp\n",
    "from snowflake.snowpark import functions as F\n",
    "from snowflake.snowpark import types as T\n",
    "\n",
    "import pandas as pd\n",
    "from sklearn.metrics import mean_squared_error\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "from statsmodels.graphics.tsaplots import plot_acf\n",
    "from statsmodels.tsa.arima.model import ARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trips_table_name = state_dict['trips_table_name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snowdf = session.table(trips_table_name)\n",
    "\n",
    "top_stations = snowdf.filter(F.col('START_STATION_ID').is_not_null()) \\\n",
    "                    .groupBy('START_STATION_ID') \\\n",
    "                    .count() \\\n",
    "                    .sort('COUNT', ascending=False) \\\n",
    "                    .toPandas()['START_STATION_ID'].values.tolist()\n",
    "\n",
    "df = snowdf.filter(F.col('START_STATION_ID') == top_stations[0]) \\\n",
    "      .withColumn('DATE', \n",
    "                  F.call_builtin('DATE_TRUNC', ('DAY', F.col('STARTTIME')))) \\\n",
    "      .groupBy('DATE') \\\n",
    "      .count() \\\n",
    "      .sort('DATE').toPandas()\n",
    "\n",
    "plt.figure(figsize=(15, 8))\n",
    "ax = sns.lineplot(x='DATE', y='COUNT', data=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use('seaborn-darkgrid')\n",
    "plt.rc(\"figure\", figsize=(10, 7))\n",
    "plot_acf(df['COUNT'], lags=400)\n",
    "plt.xlabel('Lags', fontsize=12)\n",
    "plt.ylabel('Autocorrelation', fontsize=12)\n",
    "plt.title('Autocorrelation of Trip Count Seasonality', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can definitely see the strong annual seasonality.  Lets look closer at the daily and weekly lag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc(\"figure\", figsize=(10, 7))\n",
    "plot_acf(df['COUNT'], lags=[1, 7, 30, 60, 90, 365])\n",
    "plt.xlabel('Lags', fontsize=12)\n",
    "plt.ylabel('Autocorrelation', fontsize=12)\n",
    "plt.title('Autocorrelation of Trip Count Seasonality', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline Model\n",
    "Lets build a baseline with ARIMA since we already have statsmodels imported"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima.model import ARIMA\n",
    "model=ARIMA(df['COUNT'],order=(1,1,1))\n",
    "history=model.fit()\n",
    "df['HISTORY']=history.predict(start=0, end=len(df))\n",
    "plt.figure(figsize=(15, 8))\n",
    "df1 = pd.melt(df, id_vars=['DATE'], value_vars=['COUNT', 'HISTORY'])\n",
    "ax = sns.lineplot(x='DATE', y='value', hue='variable', data=df1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"P-Value = \", adfuller(df['COUNT'].dropna(), autolag = 'AIC')[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a p-value greater than .05 we know that the trend is non-trivial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rolling_mean = df['COUNT'].rolling(window = 7).mean()\n",
    "df['STATIONARY'] = rolling_mean - rolling_mean.shift()\n",
    "ax1 = plt.subplot()\n",
    "df['STATIONARY'].plot(title='Differenced');\n",
    "ax2 = plt.subplot()\n",
    "df['COUNT'].plot(title='original')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"P-Value = \", adfuller(df['STATIONARY'].dropna(), autolag = 'AIC')[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By differencing with the rolling mean we can stationarize the series.  In order to account for this trend we can create an exogenous signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "# exog = pd.DataFrame(pd.PeriodIndex(df['DATE'], freq='D')).set_index('DATE')\n",
    "# exog['f1'] = np.sin(2 * np.pi * exog.index.dayofyear / 365.25)\n",
    "# exog['f2'] = np.cos(2 * np.pi * exog.index.dayofyear / 365.25)\n",
    "# exog['f3'] = np.sin(4 * np.pi * exog.index.dayofyear / 365.25)\n",
    "# exog['f4'] = np.cos(4 * np.pi * exog.index.dayofyear / 365.25)\n",
    "# exog = exog.reset_index()\n",
    "# exog = exog.drop(columns=['DATE'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...or use we can use the `trend=` flag in ARIMA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model=ARIMA(endog=df['COUNT'][:-365], trend='ct', order=(1,0,1))\n",
    "history=model.fit()\n",
    "df['HISTORY']=history.predict(start=0, end=len(df))\n",
    "plt.figure(figsize=(15, 8))\n",
    "df1 = pd.melt(df[:-365], id_vars=['DATE'], value_vars=['COUNT', 'HISTORY'])\n",
    "ax = sns.lineplot(x='DATE', y='value', hue='variable', data=df1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "dftest=pd.DataFrame()\n",
    "dftest['DATE']=df['DATE'][-365:]\n",
    "dftest['COUNT']=df['COUNT'][-365:]\n",
    "dftest['FORECAST']=history.forecast(steps=365)\n",
    "plt.figure(figsize=(15, 8))\n",
    "df1 = pd.melt(dftest, id_vars=['DATE'], value_vars=['COUNT', 'FORECAST'])\n",
    "ax = sns.lineplot(x='DATE', y='value', hue='variable', data=df1)\n",
    "\n",
    "error = mean_squared_error(dftest['COUNT'], dftest['FORECAST'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ARIMA Error is: \"+str(error))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "authors": [
   {
    "name": "cforbe"
   }
  ],
  "kernelspec": {
   "display_name": "snowpark_0110:Python",
   "language": "python",
   "name": "conda-env-snowpark_0110-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "msauthor": "trbye"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}