Demand_Forecasting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.tsa.api as smt\n",
    "\n",
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from keras.models import Sequential\n",
    "from keras.layers import LSTM, Dense, Dropout, Flatten, TimeDistributed, GRU\n",
    "from keras.layers.convolutional import Conv1D, MaxPooling1D\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "pd.options.display.max_rows = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('Air_Traffic_Passenger_Statistics.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dropping nulls that I found in columns - only a few\n",
    "df.dropna(inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting one airline's flights against time to see how the time series look\n",
    "\n",
    "df_explore_virgin = df[(df['Operating Airline'] == 'Virgin Atlantic')  & (df['Activity Type Code'] == 'Enplaned')]\n",
    "\n",
    "# Plot of passengers leaving SFO on Virgin Atlantic\n",
    "plt.figure(figsize = (15,8))\n",
    "plt.plot(df_explore_virgin['Activity Period'], df_explore_virgin['Passenger Count']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_explore = df[df['Activity Type Code'] == 'Enplaned'].groupby('Activity Period').sum()\n",
    "\n",
    "df_explore = df_explore[df_explore.index >= 201800]\n",
    "\n",
    "# Zoomed-in version shows that it is interpreting the difference between \"201812\" and \"201901\" incorrectly\n",
    "plt.figure(figsize = (15,8))\n",
    "plt.plot(df_explore.index, df_explore['Passenger Count'])\n",
    "plt.savefig('img/data_cleaning_dates.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['date'] = pd.to_datetime(df['Activity Period'], format = '%Y%m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Eliminating rows that account for passengers in layovers. Don't care about those\n",
    "df.set_index('Activity Type Code', inplace = True)\n",
    "df.drop(index = 'Thru / Transit', inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Changing the index to dates so that I can use time-series unique functions like pd.resample()\n",
    "df.reset_index(inplace = True)\n",
    "\n",
    "# Dropping the no longer useful column that we derived dates from\n",
    "df.drop(columns = 'Activity Period', inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cleaning up the name of the 'United Airlines' that is pre-2013\n",
    "df.replace('United Airlines - Pre 07/01/2013', 'United Airlines', inplace = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next up: I need to group by the proper types so I can get \"enplaned\" and \"deplaned\" trends, and conduct separate time series analysis on them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_explore_virgin = df[(df['Operating Airline'] == 'Virgin Atlantic')  & (df['Activity Type Code'] == 'Enplaned')]\n",
    "\n",
    "# Improved plot of Virgin Atlantic after datetime was fixed\n",
    "plt.figure(figsize = (12,6))\n",
    "plt.plot(df_explore_virgin.index, df_explore_virgin['Passenger Count']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enplaned = df[df['Activity Type Code'] == 'Enplaned'].groupby('date').sum()\n",
    "\n",
    "# Improved plot of total air traffic after datetime was fixed\n",
    "\n",
    "plt.figure(figsize = (15,5))\n",
    "sns.lineplot(enplaned.index, enplaned['Passenger Count'] / 1000000)\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.savefig('img/sfo_traffic.png', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (15,8))\n",
    "sns.lineplot(enplaned.index, enplaned['Passenger Count'] / 1000000)\n",
    "plt.savefig('img/cleaned_time_series.png', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operating Airline vs Published Airline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A company like SkyWest Airline is never a Published Airline - only an Operating Airline\n",
    "df.loc[df['Operating Airline'] == 'SkyWest Airlines', 'Published Airline'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Operating Airline'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Published Airline'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shows that Operating Airline is different from Published Airline ~2000 times: 1/10 of the time in my df\n",
    "df[df['Operating Airline'] != df['Published Airline']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**I will use Operating Airlines since that indicates who is actually flying the flight. The Published Airline is the one who is marketing the flight, and Operating is the one who is actually flying it. Individuals know that the flight is being operated by SkyWest Airlines, for instance, even if they bought the ticket from United.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Autocorrelation across all airlines is very close to 1\n",
    "enplaned['Passenger Count'].autocorr(lag = 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We definitely have autocorrelation - peaking at the 12 month mark and the 1 month mark\n",
    "plt.figure(figsize = (15,6))\n",
    "pd.plotting.autocorrelation_plot(enplaned['Passenger Count'])\n",
    "plt.xlim(0,60);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check for stationarity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aiming for an autocorrelation that is approximately 0, and not more negative that -0.5, per https://people.duke.edu/~rnau/411arim2.htm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(enplaned['Passenger Count'].autocorr())\n",
    "print(enplaned['Passenger Count'].diff().autocorr()) # Autocorrelation approx == 0\n",
    "print(enplaned['Passenger Count'].diff().diff().autocorr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(enplaned.std())\n",
    "print(enplaned.diff().std()) # Std is at its lowest\n",
    "print(enplaned.diff().diff().std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both of the above tests confirm that only one layer of differencing is needed in our model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARIMA guides for further work, if needed:\n",
    "\n",
    "Textbooks:\n",
    "\n",
    "https://people.duke.edu/~rnau/seasarim.htm\n",
    "\n",
    "http://www.statsoft.com/Textbook/Time-Series-Analysis#1general\n",
    "\n",
    "Articles:\n",
    "\n",
    "https://machinelearningmastery.com/sarima-for-time-series-forecasting-in-python/\n",
    "\n",
    "https://medium.com/@josemarcialportilla/using-python-and-auto-arima-to-forecast-seasonal-time-series-90877adff03c\n",
    "\n",
    "Documentation examples:\n",
    "\n",
    "https://www.statsmodels.org/stable/examples/notebooks/generated/statespace_sarimax_stata.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_enplaned = df[df['Activity Type Code'] == 'Enplaned']\n",
    "df_deplaned = df[df['Activity Type Code'] == 'Deplaned']\n",
    "\n",
    "df_enplaned = df_enplaned.groupby(by = ['date', 'Operating Airline']).sum()\n",
    "df_deplaned = df_deplaned.groupby(by = ['date', 'Operating Airline']).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_deplaned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_grouped = df_enplaned.merge(df_deplaned,\n",
    "                               left_index = True,\n",
    "                               right_index = True,\n",
    "                               how = 'left',\n",
    "                               suffixes = (' Enplaned', ' Deplaned'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Why the small difference between the merged and original values?\n",
    "# Difference is so small, it's probably not worth investigating. Couldn't find a 1550 passenger row\n",
    "print(df_grouped['Passenger Count Enplaned'].sum(), df_enplaned['Passenger Count'].sum())\n",
    "print(df_grouped['Passenger Count Deplaned'].sum(), df_deplaned['Passenger Count'].sum())\n",
    "print(df_grouped['Passenger Count Deplaned'].sum() - df_deplaned['Passenger Count'].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_grouped['Passenger Count Difference'] = df_grouped['Passenger Count Deplaned'] - \\\n",
    "                                           df_grouped['Passenger Count Enplaned']\n",
    "\n",
    "df_grouped['Passenger Count %'] = df_grouped['Passenger Count Deplaned'] / \\\n",
    "                                  df_grouped['Passenger Count Enplaned'] * 100\n",
    "\n",
    "df_grouped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_grouped.reset_index(inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_enplaned = df_grouped.groupby('Operating Airline').sum()['Passenger Count Enplaned'].sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cutoff seems to be around 250,000 total passengers, for major airlines\n",
    "total_enplaned[:-1].plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can see the most important airways\n",
    "list(zip(total_enplaned.index, list(total_enplaned)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finds the names of the top 15 airlines by quantity of passengers leaving SFO\n",
    "\n",
    "top15_list = list(df_grouped.groupby('Operating Airline').\n",
    "                  sum().\n",
    "                  sort_values('Passenger Count Enplaned', ascending = False)[:15].\n",
    "                  index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top15_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time Series (full airport's departures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_airport = df_deplaned.groupby('date').sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train/Test Split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_airport.diff().plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taking the log of the train data for training the model on multiplicative data\n",
    "train = dp_airport.loc[dp_airport.index.year <= 2014, 'Passenger Count']\n",
    "test = dp_airport.loc[dp_airport.index.year >= 2015, 'Passenger Count']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Visualization, pre-SARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1 diff seems to work really well for this data\n",
    "dp_airport.diff()['Passenger Count'].autocorr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define helper plot function for visualization\n",
    "# Needs to be fixed so that data.plot either works (damaged by FB Prophet import)\n",
    "def make_plots(data, lags=None):\n",
    "    layout = (1, 3)\n",
    "    raw  = plt.subplot2grid(layout, (0, 0))\n",
    "    acf  = plt.subplot2grid(layout, (0, 1))\n",
    "    pacf = plt.subplot2grid(layout, (0, 2))\n",
    "    \n",
    "    #plt.plot(data,ax=raw,figsize=(12,6))\n",
    "    data.plot(ax=raw, figsize=(12, 6))\n",
    "    smt.graphics.plot_acf(data, lags=lags, ax=acf)\n",
    "    smt.graphics.plot_pacf(data, lags=lags, ax=pacf)\n",
    "    sns.despine()\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ACF and PACF make_plots look the same from all different starting point months\n",
    "plt.figure(figsize = (15,8))\n",
    "make_plots(dp_airport.diff()[1:], lags=24);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build SARIMA pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, I'm going to have only a train/test split, with a large test set (let's say 4 years, so I can predict SFO's capacity needs 4 years into the future).\n",
    "\n",
    "Instead of testing my model's performance on predicting validation data, I'm going to pick the model with the lowest AIC (chooses the best model based on a combination of how well the model fits the data (with log-likelihood), penalized for model complexity via its degrees of freedom).\n",
    "\n",
    "Then, once I have the model with the lowest AIC, I'm going to use that model to predict my test set. (https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAXResults.predict.html)\n",
    "\n",
    "Once I have my predictions, I'll calculate the residuals of the predictions, divided by the true values. **(Note: I should only be scoring the performance of my model on predictions it makes starting 2 years after it sees data, since the ideal model deals with 2 previous periods on a 12 month lag)** This will give me a good measure of how much error to expect in my capacity planning (Ex: 5%). Useful for presentation purposes on what to expect from the model.\n",
    "\n",
    "Then I will retrain my model with the most recent data, and attempt to predict the values for 2019-23. The expected error is the same as what we had from my residuals in my 2015-19 prediction."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Grid Searching my way to the ideal SARIMA model\n",
    "sar_aic = {}\n",
    "\n",
    "for p in range(3):\n",
    "    for d in range(2):\n",
    "        for q in range(3):\n",
    "            for P in range(3):\n",
    "                for D in range(2):\n",
    "                    for Q in range(3):\n",
    "                        for trend in ['n','c','t','ct']:\n",
    "                            sar = SARIMAX(train, order=(p,d,q), seasonal_order=(P,D,Q,12), trend=trend,\n",
    "                                          enforce_stationarity=False, enforce_invertibility=False).fit()\n",
    "\n",
    "                            sar_aic[f\"SARIMA({p},{d},{q})({P},{D},{Q},12), {t} trend\"] = sar.aic\n",
    "                    "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "min(sar_aic, key = sar_aic.get)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sar = SARIMAX(train, order=(1,0,2), seasonal_order=(0,1,2,12), trend='t',\n",
    "              enforce_stationarity=False, enforce_invertibility=False).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "make_plots(sar.resid, lags = 48)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### The above plot shows that we can still get predictive power out of the first order ACF and PACF\n",
    "\n",
    "### Below - a new model defined with a second order of differencing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note: This cell had trend of \"t\" instead of \"n\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_sarima = SARIMAX(train, order=(1,1,2), seasonal_order=(0,1,2,12), trend='t',\n",
    "                      enforce_stationarity=False, enforce_invertibility=False).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should show ACF and PACF plot - look at make_plots to fix broken function due to FB Prophet\n",
    "make_plots(best_sarima.resid, lags = 48)\n",
    "plt.savefig('img/acf_and_pacf.png', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "print(min(sar_aic, key = sar_aic.get), sar_aic[min(sar_aic, key = sar_aic.get)])\n",
    "print('SARIMA(1,1,2)(0,1,2,12), t trend', best_sarima.aic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model predictions start 2 years out, since model incorporates 2 MA elements of 12 month seasons\n",
    "#plt.figure(figsize=(15,8))\n",
    "train[25:].plot()\n",
    "best_sarima.predict(start = 25).plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,5))\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.plot(dp_airport[25:] / 1000000)\n",
    "plt.plot(best_sarima.predict(start=25, end=len(train)+len(test)) / 1000000)\n",
    "plt.savefig('img/sarima_fit.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Starts prediction at the end of train data, through the full dataset (train + test)\n",
    "test_prediction = best_sarima.predict(start=len(train), end=(len(train)+len(test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.figure(figsize=(12,7))\n",
    "test.plot()\n",
    "test_prediction.plot()\n",
    "plt.savefig('img/test_prediction.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arima_residuals = test_prediction - test\n",
    "plt.scatter(test_prediction.index, arima_residuals);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mae_residuals(start, end, test):\n",
    "    '''\n",
    "    Finds residuals based on \"best_sarima\" model, and the end point defined relative to the training data\n",
    "    '''\n",
    "    test_predictions = best_sarima.predict(start=len(train)+start, end=len(train)+end)\n",
    "    residuals = test_predictions - test\n",
    "    answer = abs(residuals).sum() / test[start:end].sum() * 100 # Calculate residual MAE for the proper slice of time\n",
    "    return answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# My predictions can be overall expected to be ~5.25% off from the actual values\n",
    "abs(arima_residuals).sum() / test.sum() * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***My ARIMA predictions can be expected to be ~5.25% off from the actual values - projection below***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAE by year, to show increasing error over time for model's predictions\n",
    "print(mae_residuals(0,12, test))\n",
    "print(mae_residuals(12, 24, test))\n",
    "print(mae_residuals(24, 36, test))\n",
    "print(mae_residuals(36, 48, test))\n",
    "\n",
    "print(mae_residuals(0,48, test)) # Overall residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***The above are the errors that I can expect for every year of subsequent predictions***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_sarima.plot_diagnostics(figsize = (15,8));\n",
    "plt.savefig('img/sarima_diagnostics.png', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_sarima = SARIMAX(train, order=(1,1,2), seasonal_order=(0,1,2,12), trend='t',\n",
    "                       enforce_stationarity=False, enforce_invertibility=False).fit()\n",
    "\n",
    "test_sarima_prediction = test_sarima.predict(start=len(train), end=len(train)+96)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_airport.plot()\n",
    "test_sarima_prediction.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final SARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: This final_sarima is tweaked from the ideal one above, because the ARIMA was projecting a declining\n",
    "# SFO attendance, probably due to the two most recent years of data that it erroneously extrapolates from\n",
    "final_sarima = SARIMAX(dp_airport, order=(1,1,2), seasonal_order=(1,1,2,12), trend='n',\n",
    "                       enforce_stationarity=False, enforce_invertibility=False).fit()\n",
    "\n",
    "final_sarima_prediction = final_sarima.predict(start=len(dp_airport), end=len(dp_airport)+48)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(13,7))\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.plot(dp_airport.index, dp_airport / 1000000)\n",
    "plt.plot(final_sarima_prediction.index, final_sarima_prediction / 1000000)\n",
    "plt.savefig('img/sarima_projection.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the final diagnostic plot of my best SARIMA\n",
    "final_sarima.plot_diagnostics(figsize = (15,8))\n",
    "plt.savefig('img/sarima_assumptions.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(final_sarima_prediction) -  np.max(dp_airport)) # Difference in \n",
    "print(np.max(final_sarima_prediction) / np.max(dp_airport))\n",
    "\n",
    "np.set_printoptions(suppress=True)\n",
    "print(final_sarima_prediction[-12:].sum() / dp_airport[-12:].sum()) # Percentage traffic growth after 4 years\n",
    "print((final_sarima_prediction[-12:].sum() - dp_airport[-12:].sum()) / 4) # Avg traffic growth each year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**My SARIMA model predicts almost no growth in passenger count - only 2% growth in 4 years. My SARIMA model takes the last 2 year's leveling off of the trend very seriously**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Facebook Prophet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Room for improvement - include cross-validation over time with FBProphet**\n",
    "https://facebook.github.io/prophet/docs/diagnostics.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fbprophet import Prophet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting up data in proper format for Prophet - and dividing into millions\n",
    "fb_train = (train / 1000000).reset_index()\n",
    "fb_train.columns = ['ds','y']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fb_model = Prophet(seasonality_mode='multiplicative')\n",
    "fb_model.fit(fb_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb_test = test.reset_index()\n",
    "fb_test.columns = ['ds','y']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecast = fb_model.predict(fb_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot of FBProphet performance (multiplicative) vs the test set\n",
    "fb_model.plot(forecast)\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "#plt.set_yticklabels(labels = [1.25,1.5,1.75,2.0,2.25,2.5,2.75,3.0])\n",
    "plt.grid(b=True)\n",
    "plt.gca().get_lines()[1].set_color(\"orange\")\n",
    "plt.plot(test.index, test / 1000000)\n",
    "plt.savefig('img/fb_fit.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test / 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb_pred = forecast['yhat'] # Separating the prediction\n",
    "fb_pred.index = test.index # Allowing them to be subtracted with the same index\n",
    "fb_mae = np.mean(abs(fb_pred - (test / 1000000)) / (test / 1000000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb_mae"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**FB Prophet's model is only *slightly* better than SARIMA - at 5.4% MAE rather than SARIMA's 5.5%**\n",
    "\n",
    "### Interpretation:\n",
    "\n",
    "Since my data is a fairly straightforward time series without much outside data, it is performing in SARIMA's strong point, without using the advantages of more sophisticated models. I can expect that the LSTM will have a similar performance, since it may be an actual anomaly in growth at SFO that contributed to the 5.5% MAE from prediction to actual values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb_model.plot_components(forecast);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb_dp_airport = (dp_airport / 1000000).reset_index()\n",
    "fb_dp_airport.columns = ['ds','y']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Final FB Prophet model to predict 4 years into the future\n",
    "fb_model_final = Prophet(seasonality_mode='multiplicative').fit(fb_dp_airport)\n",
    "\n",
    "future = fb_model_final.make_future_dataframe(periods=48, freq='M') # Making new rows for future\n",
    "future[-48:] += pd.DateOffset(1) # +1 on dates to change from 3/31 to 4/01, for example. Uniform w/ previous data\n",
    "\n",
    "fb_future_pred = fb_model_final.predict(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting predictions into the future\n",
    "fb_model_final.plot(fb_future_pred[-49:])\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.gca().get_lines()[1].set_color(\"orange\")\n",
    "plt.savefig('img/fb_projection.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predicted needed max capacity of SFO in the year 2023\n",
    "print(np.max(fb_future_pred['yhat']))\n",
    "print(np.max(test))\n",
    "\n",
    "print(fb_future_pred['yhat'][-12:].sum() / dp_airport[-12:].sum()) # Percentage traffic growth after 4 years\n",
    "print((fb_future_pred['yhat'][-12:].sum() - dp_airport[-12:].sum()) / 4) # Avg traffic growth each year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LSTM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking to exceed the performance of my SARIMA and FB Prophet models in predicting the test set, after training on the train set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Starting the training data at the beginning of the year\n",
    "lstm_train = train[6:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data preprocessing in order to get rows of 12 months as x, and the result as y\n",
    "# Credit to Jason Brownlee: https://machinelearningmastery.com/how-to-develop-lstm-models-for-time-series-forecasting/\n",
    "\n",
    "def split_sequence(sequence, n_steps):\n",
    "    X, y = list(), list()\n",
    "    for i in range(len(sequence)):\n",
    "        # find the end of this pattern\n",
    "        end_ix = i + n_steps\n",
    "        # check if we are beyond the sequence\n",
    "        if end_ix > len(sequence)-1:\n",
    "            break\n",
    "        # gather input and output parts of the pattern\n",
    "        seq_x, seq_y = sequence[i:end_ix], sequence[end_ix]\n",
    "        X.append(seq_x)\n",
    "        y.append(seq_y)\n",
    "    return np.array(X), np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = split_sequence(lstm_train, 12)\n",
    "\n",
    "X = X.reshape((96,12,1)) # Adding an artificial third dimension to fit into my LSTM model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lstm = Sequential()\n",
    "\n",
    "# Start with simple model - 1 or 2 layers. Fiddle with dropouts, batch size, epochs\n",
    "# Length of training data, how far back the window goes, # of features to predict\n",
    "lstm.add(LSTM(64, activation='relu', input_shape=(12,1)))\n",
    "#lstm.add(LSTM(32, activation='relu', input_shape=(32,1)))\n",
    "lstm.add(Dense(10))\n",
    "lstm.add(Dropout(0.3))\n",
    "#lstm.add(Dense(10))\n",
    "#lstm.add(Dropout(0.3))\n",
    "lstm.add(Dense(1))\n",
    "\n",
    "\n",
    "lstm.compile(optimizer = 'adam', loss = 'mae', )\n",
    "# lstm.fit(X, y, shuffle=False, batch_size=8, epochs=20)\n",
    "lstm.fit(X, y, validation_data=(X[-24:],y[-24:]), shuffle=False, batch_size=32, epochs=100)\n",
    "# Batch size of 8 or 16 - increments of 2^x\n",
    "\n",
    "# Specify batch size when I fit model - larger is more accurate but slower"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lstm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNN - LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Testing out CNN-LSTMs\n",
    "\n",
    "n_features = 1\n",
    "# Splitting it into seasons - to see if there are patterns\n",
    "n_seq = 4\n",
    "# Considering 1 year's worth of data. Should experiment with 24 months - based on SARIMA\n",
    "n_steps = 3\n",
    "X = X.reshape((X.shape[0], n_seq, n_steps, n_features))\n",
    "\n",
    "# define model\n",
    "model = Sequential()\n",
    "model.add(TimeDistributed(Conv1D(filters=64, kernel_size=1, activation='relu'), input_shape=(None, n_steps, n_features)))\n",
    "model.add(TimeDistributed(MaxPooling1D(pool_size=2)))\n",
    "model.add(TimeDistributed(Flatten()))\n",
    "model.add(LSTM(64, activation='relu'))\n",
    "model.add(Dense(1))\n",
    "model.compile(optimizer='adam', loss='mae')\n",
    "# fit model\n",
    "model.fit(X, y, validation_data=(X[-24:],y[-24:]), shuffle=False, batch_size=16, epochs=100)\n",
    "# demonstrate prediction\n",
    "# x_input = array([60, 70, 80, 90])\n",
    "# x_input = x_input.reshape((1, n_seq, n_steps, n_features))\n",
    "# yhat = model.predict(x_input, verbose=0)\n",
    "# print(yhat)\n",
    "\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "cnn_lstm_test = np.array(test).reshape(1, n_seq, n_steps, n_features)\n",
    "model.predict(cnn_lstm_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GRU (trained on train data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = split_sequence(lstm_train, 12)\n",
    "X = X.reshape(X.shape[0], X.shape[1], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru = Sequential()\n",
    "\n",
    "# Start with simple model - 1 or 2 layers. Fiddle with dropouts, batch size, epochs\n",
    "# Length of training data, how far back the window goes, # of features to predict\n",
    "gru.add(GRU(128, activation='relu', input_shape=(12,1), return_sequences=True))\n",
    "gru.add(GRU(128, activation='relu', input_shape=(32,1)))\n",
    "gru.add(Dense(64))\n",
    "# gru.add(Dropout(0.08))\n",
    "gru.add(Dense(32))\n",
    "gru.add(Dense(16))\n",
    "# gru.add(Dropout(0.08))\n",
    "gru.add(Dense(1))\n",
    "\n",
    "\n",
    "gru.compile(optimizer = 'adam', loss = 'mae', )\n",
    "# lstm.fit(X, y, shuffle=False, batch_size=8, epochs=20)\n",
    "gru.fit(X, y, validation_data=(X[-24:],y[-24:]), shuffle=False, batch_size=32, epochs=500)\n",
    "# Batch size of 8 or 16 - increments of 2^x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_X, test_y = split_sequence(test, 12)\n",
    "test_X = test_X.reshape(test_X.shape[0], test_X.shape[1], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shows how well the model fits the train data\n",
    "gru_train_pred = gru.predict(X)\n",
    "\n",
    "plt.plot(y)\n",
    "plt.plot(gru_train_pred);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  vvv Old prediction that may be wrong based on my analysis"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Predicts the end of the train data, into the beginning of test data (fills gap from test prediction)\n",
    "gru_start_test_pred = gru.predict(X[-12:])\n",
    "\n",
    "# Shows how well the model fits the test data\n",
    "gru_test_pred = gru.predict(test_X)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "print(gru_start_test_pred.shape, gru_test_pred.shape)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "gru_test_pred = np.append(gru_start_test_pred, gru_test_pred, axis=0)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "gru_test_pred"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Assigning dates to gru_test_pred to allow it to be plotted\n",
    "gru_test_pred = pd.DataFrame(gru_test_pred)\n",
    "gru_test_pred.index = dp_airport[-39:].index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ^^^ Old prediction that may be wrong based on my analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Feeding the last 12 months of test to my GRU now that it's been trained. My model will predict 4 years' future\n",
    "\n",
    "def rolling_predict(starting_data, keras_model, window_size, n_months_to_predict):\n",
    "    '''\n",
    "    starting_data = Starting x_month length seed to predict the next value (Ex: 12 vals to predict 13th)\n",
    "    \n",
    "    Make a rolling window prediction such that my neural nets can predict future times off their past predictions\n",
    "    '''\n",
    "    \n",
    "    # Seed the start of the rolling window\n",
    "    rolling_window = starting_data[-window_size:].reshape(1,window_size,1)\n",
    "    \n",
    "    model_predictions = []\n",
    "    \n",
    "    for i in range(n_months_to_predict): # Predicting 4 years\n",
    "        # Predicting off of the last n number of predictions according to window size, indexing the array\n",
    "        next_pred = keras_model.predict(rolling_window[:,-window_size:,:])\n",
    "        # Pulling the next_pred out of its 1x1 matrix to put into its own list\n",
    "        model_predictions.append(next_pred[0,0])\n",
    "        # Moving rolling window in the proper dimension of the 3D array needed for keras (2nd dimension)\n",
    "        rolling_window = np.append(rolling_window, next_pred.reshape(1,1,1), axis=1)\n",
    "        \n",
    "    return model_predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru_test_predictions = rolling_predict(y[-12:], gru, 12, 51)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assigning dates to gru_test_pred to allow it to be plotted\n",
    "gru_test_predictions = pd.DataFrame(gru_test_predictions)\n",
    "gru_test_predictions.index = dp_airport[-51:].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(test)\n",
    "plt.plot(gru_test_predictions);"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "train_mae = []\n",
    "\n",
    "# Iterate through to calculate the individual predictions vs actual values\n",
    "for i in range(len(test_y)):\n",
    "    mae.append(abs(gru_train_pred[i] - y[i]) / y[i])\n",
    "    \n",
    "np.mean(mae)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "mae = []\n",
    "\n",
    "# Iterate through to calculate the individual predictions vs actual values\n",
    "for i in range(len(test_y)):\n",
    "    mae.append(abs(gru_test_predictions[i] - test_y[i]) / test_y[i])\n",
    "    \n",
    "np.mean(mae)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAE is 5.7%\n",
    "np.mean(abs(gru_test_predictions.iloc[:,0] - test) / test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(13,7))\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.plot(dp_airport.index, dp_airport['Passenger Count'] / 1000000)\n",
    "plt.plot(gru_test_predictions.index, gru_test_predictions / 1000000)\n",
    "plt.savefig('img/gru_test_prediction.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru_predictions = rolling_predict(test_y, gru, 12, 48)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting it in a format that is able to be graphed\n",
    "future = future[-48:]\n",
    "future['Predictions'] = gru_predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graphing my GRU neural net\n",
    "# plt.figure(figsize=(13,7))\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.plot(dp_airport.index, dp_airport / 1000000)\n",
    "plt.plot(future['ds'], future['Predictions'] / 1000000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(gru_predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GRU 2 (trained on dp_airport)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refitting the GRU on all the data from dp_airport, so that my neural net follows the same model as the previous models (train, use test to estimate MAE, and then retrain on all data to predict final values). Unclear whether this model is best at predicting, since it gives the most extreme capacity prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru2_X, gru2_y = split_sequence(dp_airport['Passenger Count'], 12)\n",
    "gru2_X = gru2_X.reshape(gru2_X.shape[0], gru2_X.shape[1], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru2 = Sequential()\n",
    "\n",
    "# Start with simple model - 1 or 2 layers. Fiddle with dropouts, batch size, epochs\n",
    "# Length of training data, how far back the window goes, # of features to predict\n",
    "gru2.add(GRU(128, activation='relu', input_shape=(12,1), return_sequences=True))\n",
    "gru2.add(GRU(128, activation='relu', input_shape=(32,1)))\n",
    "gru2.add(Dense(64))\n",
    "# gru2.add(Dropout(0.08))\n",
    "gru2.add(Dense(32))\n",
    "gru2.add(Dense(16))\n",
    "# gru2.add(Dropout(0.08))\n",
    "gru2.add(Dense(1))\n",
    "\n",
    "\n",
    "gru2.compile(optimizer = 'adam', loss = 'mae', )\n",
    "# lstm.fit(X, y, shuffle=False, batch_size=8, epochs=20)\n",
    "gru2.fit(gru2_X, gru2_y, validation_data=(gru2_X[-24:],gru2_y[-24:]), shuffle=False, batch_size=16, epochs=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru2_train_pred = gru.predict(gru2_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(gru2_y)\n",
    "plt.plot(gru2_train_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru2_predictions = rolling_predict(gru2_y, gru2, 12, 48)\n",
    "\n",
    "future['Predictions_2'] = gru2_predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gru2 was retrained on all available time series data (minus validation), while gru was trained only on train set\n",
    "# plt.figure(figsize=(13,7))\n",
    "# plt.title('Departures (in millions)', fontsize=30)\n",
    "# plt.xlabel('')\n",
    "# plt.ylabel('')\n",
    "# plt.xticks(fontsize=17)\n",
    "# plt.yticks(fontsize=17)\n",
    "# plt.grid(b=True)\n",
    "plt.plot(dp_airport.index, dp_airport / 1000000)\n",
    "plt.plot(future['ds'], future['Predictions_2'] / 1000000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GRU 3 (no validation data separate on the test data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru3_X, gru3_y = split_sequence(dp_airport['Passenger Count'], 12)\n",
    "gru3_X = gru3_X.reshape(gru3_X.shape[0], gru3_X.shape[1], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gru3 = Sequential()\n",
    "\n",
    "# Start with simple model - 1 or 2 layers. Fiddle with dropouts, batch size, epochs\n",
    "# Length of training data, how far back the window goes, # of features to predict\n",
    "gru3.add(GRU(128, activation='relu', input_shape=(12,1), return_sequences=True))\n",
    "gru3.add(GRU(128, activation='relu', input_shape=(32,1)))\n",
    "gru3.add(Dense(64))\n",
    "# gru3.add(Dropout(0.08))\n",
    "gru3.add(Dense(32))\n",
    "gru3.add(Dense(16))\n",
    "# gru3.add(Dropout(0.08))\n",
    "gru3.add(Dense(1))\n",
    "\n",
    "\n",
    "gru3.compile(optimizer = 'adam', loss = 'mae', )\n",
    "# lstm.fit(X, y, shuffle=False, batch_size=8, epochs=20)\n",
    "gru3.fit(gru3_X, gru3_y, shuffle=False, batch_size=32, epochs=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gru3_predictions = rolling_predict(gru3_y, gru3, 12, 48)\n",
    "\n",
    "future['Predictions_3'] = gru3_predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(13,7))\n",
    "plt.title('Departures (in millions)', fontsize=30)\n",
    "plt.xlabel('')\n",
    "plt.ylabel('')\n",
    "plt.xticks(fontsize=17)\n",
    "plt.yticks(fontsize=17)\n",
    "plt.grid(b=True)\n",
    "plt.plot(dp_airport.index, dp_airport / 1000000)\n",
    "plt.plot(future['ds'], future['Predictions_3'] / 1000000)\n",
    "plt.savefig('img/gru_future_prediction.png', bbox_inches = 'tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The first GRU is for tweaking the model hyperparameters and architecture, and getting the test score. GRU2 and GRU3 are both for predicting the future, but GRU2 withholds the last 2 years as validation, while GRU3 is trained on everything"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Growth rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Printing out the percentage of total growth from 4 years out to \n",
    "print(future.iloc[-12:, 1].sum() / dp_airport[-12:].sum()) # GRU\n",
    "print(future.iloc[-12:, 2].sum() / dp_airport[-12:].sum()) # GRU2\n",
    "print(future.iloc[-12:, 3].sum() / dp_airport[-12:].sum()) # GRU3\n",
    "\n",
    "\n",
    "print((future.iloc[-12:, 1].sum() - dp_airport[-12:].sum()) / 4)\n",
    "print((future.iloc[-12:, 2].sum() - dp_airport[-12:].sum()) / 4)\n",
    "print((future.iloc[-12:, 3].sum() - dp_airport[-12:].sum()) / 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final GRU decision:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I made three GRU models to see which one would predict the future the best. The fact that I need to retrain a new neural net on all the data makes it more tricky than a simpler modeling method like SARIMA or FB Prophet. In the end, I normally would have chosen to use GRU2, which is the GRU that retrains on all the data, but holds out the last 2 years as validation data that the model will be scored against, just like the original model that was trained.\n",
    "\n",
    "Otherwise, if I attempt to train on all the data through 2019, without holding out a validation set, I will be training a different neural net, and it would:\n",
    "\n",
    "1. Possibly overfit\n",
    "2. More importantly: *Be a different model than the one I had trained on, and gotten a MAE for*\n",
    "\n",
    "Though holding 2017-19 as a validation set means my GRU has less data than SARIMA and FB Prophet did when they retrained, it is necessary for the integrity of my neural net."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### However:\n",
    "After looking at the results of my GRUs, you can see that the dataset where the final 2 years are held out predict *dramatically* higher growth. This makes me think that the last 2 years are important in order to factor in the tapering-out of the airport traffic growth trend.\n",
    "\n",
    "After looking at the smoothness of the graph, I could be confident that the model was not overfitting, since it was not even capturing the mini-seasonality that past years consistently show. This made me more comfortable choosing GRU3 (the GRU that was retrained on all data from '05-'19), over GRU2 (the GRU that was retrained on '05-'17, with the last 24 months withheld as validation set for model scoring purposes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Future Work: Time series on Passenger Difference (for each of 15 airlines - to predict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I am predicting the seasonality of the difference in passengers coming in and out of SFO (not the seasonality of departures or arrivals themselves).\n",
    "\n",
    "This prediction could be useful for the airlines themselves (sales forecasting) or for consumers (general ideas of which airlines are the least busy - if you want cheap tickets or to fly on planes that have less passengers, and more overhead baggage space)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (15,8))\n",
    "plt.plot(df_grouped.loc[df_grouped['Operating Airline'] == 'United Airlines', 'date'],\n",
    "         df_grouped.loc[df_grouped['Operating Airline'] == 'United Airlines', 'Passenger Count Difference']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}