<jupyter_start><jupyter_text>Timeseries Pandas started out in the financial world, so it naturally has strong support for timeseries data.We'll look at some pandas data types and methods for manipulating timeseries data.Afterwords, we'll use [statsmodels' state space framework](http://www.statsmodels.org/stable/statespace.html) to model timeseries data.<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('default')
plt.rcParams['figure.figsize'] = (12, 6)
pd.options.display.max_rows = 10<jupyter_output><empty_output><jupyter_text>Datatypes- `pd.Timestamp` (nanosecond resolution `datetime.datetime`)- `pd.Timedelta` (nanosecond resolution `datetime.timedelta`) Pandas provides highly-performant (mostly) drop-in replacements for `datetime.datetime` (`pd.Timestamp`) and `datetime.tiemedelta` (`pd.Timedelta`).These have been tailored for efficient storage in NumPy arrays.For the most part you'll be working with `DatetimeIndex`es or `TimedeltaIndex`es, or Series / DataFrames containing these.The biggest limitation is that pandas stores `Timestamp`s at nanosecond resolution. Since they're backed by NumPy's 64-bit integer, the minimum and maximum values are<jupyter_code>pd.Timestamp.min, pd.Timestamp.max<jupyter_output><empty_output><jupyter_text>If this is a problem, [there are workarounds](http://pandas.pydata.org/pandas-docs/stable/timeseries.htmlrepresenting-out-of-bounds-spans). We'll go back to the BTS data set on flights.This time I've provided the number of flights per hour for two airports in Chicago: Midway (MDW) and O'Hare (ORD). The data go back to January 1st, 2000.<jupyter_code>df = pd.read_csv("data/flights-ts.csv.gz", index_col=0, parse_dates=True)
df.head()<jupyter_output><empty_output><jupyter_text>ResamplingResampling is similar to a groupby, but specialized for datetimes.Instead of specifying a column of values to group by, you  specify a `rule`: the desired output frequency.The original data is binned into each group created by your rule.<jupyter_code>resampler = df.resample("MS")  # MS=Month Start
resampler<jupyter_output><empty_output><jupyter_text>There's an extensive list of frequency codes: http://pandas.pydata.org/pandas-docs/stable/timeseries.htmloffset-aliases.If you examine the raw data in `df`, you'll notice that it's not at a fixed frequency.Hours where there weren't any flights just simply aren't present.This isn't a problem though; resample is perfect for going from "ragged" timeseries data to fixed-frequency data.Just like with `.groupby`, `.resample` returns a deferred object that hasn't really done any work yet.It has methods for aggregation, transformation, and general function application.<jupyter_code>resampler.sum()
resampler.sum().plot();<jupyter_output><empty_output><jupyter_text>Exercise: ResamplePlot the standard deviation for the number of flights from `MDW` and `ORD` at a weekly frequency<jupyter_code># Your solution

%load solutions/timeseries_resample.py<jupyter_output><empty_output><jupyter_text>Exercise: Resample-AggCompute the the total number of flights (sum), mean, and median flights *per Quarter*.<jupyter_code>%load solutions/timeseries_resample_agg.py<jupyter_output><empty_output><jupyter_text>Rolling, ExpandingApplying functions to windows, moving through your data. These are very similar to groupby and resample. Let's get the daily number of flights with a `resample` quick.<jupyter_code>daily = df.resample('D').sum()
daily<jupyter_output><empty_output><jupyter_text>Suppose you wanted a 30-day moving (or rolling) average.This is possible with the `.rolling` method. Like `groupby` and `resample`, this object is just going to store the information to know what subset of data to operate on next; it doesn't actually do any work yet:<jupyter_code>daily.rolling(30, center=True)<jupyter_output><empty_output><jupyter_text>The first argument is the window size.Since `daily` is at daily frequency, 30 means a 30-day window.`center=True` says to label each window with the middle-most point.To actually do work, you call a method like `.mean`;<jupyter_code>fig, ax = plt.subplots()
daily.rolling(30).mean().rename(columns=lambda x: x + " (30D MA)").plot(ax=ax, alpha=.25,
                                                                        color=['C0', 'C1'])
daily.plot(ax=ax, alpha=.25, color=['C0', 'C1'], legend=False);<jupyter_output><empty_output><jupyter_text>It's common to combine resampling and rolling.<jupyter_code>df.resample("D").sum().rolling(30).corr(pairwise=True).xs("MDW", level=1)['ORD'].plot(
    title="O'Hare : Midway cross-correlation (30D MA)", figsize=(12, 4)
);<jupyter_output><empty_output><jupyter_text>Timezonespandas can store an array of datetimes with a common timezone.Right now the index for `df` is timezone naïve, but we can convert to a timezone with `tz_convert`:<jupyter_code>df.index.tzinfo  # None, timezone naïve
df.index.tz_localize("US/Central")<jupyter_output><empty_output><jupyter_text>Timezones, as usual, are annoying to deal with.We've hit a daylight savings time issue.As the error says, 2000-04-02T02:00:00 isn't actaully a valid time in US/Central.I checked the BTS website, and these timestamps are supposed to be local time, so presumably some data was recorded incorrectly.pandas is strict by default, so it we need to tell it to ignore those errors:<jupyter_code>idx = df.index.tz_localize("US/Central", ambiguous="NaT", errors='coerce')
idx
pd.isnull(idx).sum()  # 25 bad values<jupyter_output><empty_output><jupyter_text>Notice the dtype: `datetime64[ns, US/Central]`.That means nanosecond resolution in the US/Central time zone.Once you have a datetime with timezone, you can convert timezones with `tz_convert`:<jupyter_code>idx.tz_convert("UTC")<jupyter_output><empty_output><jupyter_text>OffsetsI wish the standard library `datetime` module had something like this.Let's generate some fake data with `pd.date_range`<jupyter_code>dates = pd.date_range("2016-01-01", end="2016-12-31", freq='D')
dates<jupyter_output><empty_output><jupyter_text>There are a whole bunch of offsets available in the `pd.tseries.offsets` namespace. For example, to move 3 business days into the future:<jupyter_code>dates + pd.tseries.offsets.BDay(3)<jupyter_output><empty_output><jupyter_text>Or to move to the next month end:<jupyter_code>dates + pd.tseries.offsets.MonthEnd()<jupyter_output><empty_output><jupyter_text>Timedelta MathBeing able to add columns of dates and timedeltas turns out to be quite convenient.Let's go all the way back to our first example with flight delays from New York airports.<jupyter_code>flights = pd.read_csv("data/ny-flights.csv.gz", parse_dates=['dep', 'arr'])
flights.head()<jupyter_output><empty_output><jupyter_text>Exercise: Convert TimedeltaConvert `flights.dep_delay` and `flights.arr_delay` to timedelta dtype.- Hint: recall our type conversion methods: `pd.to_*`- Make new columns in `flights` called `dep_delay_td` and `arr_delay_td`- Check the `unit` argument for the conversion method. The delay columns are in *minutes*.<jupyter_code># Your solution

%load solutions/timeseries_timedelta.py<jupyter_output><empty_output><jupyter_text>Exercise: Timedelta MathCompute the actual time the flight left, but adding the departure time `dep` and the delay `dep_delay`.<jupyter_code>%load solutions/timeseries_departure.py<jupyter_output><empty_output><jupyter_text>Modeling TimeseriesTimeseries are an interesting problem to model.If we're lucky, we have a long history of past data that we can (maybe) use to predict the future.We can exploit regularity in the timeseries (seasonal patterns, periods of high values are typically followed by another high value, etc.) to better predict the future.Statsmodels has a nice framework for fitting timeseries models and evaluating their output.<jupyter_code>import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm<jupyter_output><empty_output><jupyter_text>Let's model Monthly flights from `ORD`.<jupyter_code>y = daily.ORD.resample("MS").sum()
y.plot();<jupyter_output><empty_output><jupyter_text>That final value is odd because it's not a complete month. Let's drop it.<jupyter_code>y = daily.ORD.resample("MS").sum().iloc[:-1]
y.head()<jupyter_output><empty_output><jupyter_text>It's common to estimate the parameters on *differenced* values.That is, make a new series $y'$ where $y_t' = y_t - y_{t-1}$. Pandas makes this simple with the `.diff` method.<jupyter_code>y_prime = y.diff()
y_prime.head()<jupyter_output><empty_output><jupyter_text>We'll drop that first NaN:<jupyter_code>y_prime = y.diff().dropna()
y_prime.plot();<jupyter_output><empty_output><jupyter_text>Think back to regular linear regression: Predict some variable $y$ with some matrix $X$:$y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 ... + \beta_p X_p + \varepsilon$When modelling timeseries, past values of $y$ make for good components of $X$.We can do this with the pandas `.shift` method:<jupyter_code>y_prime.shift()<jupyter_output><empty_output><jupyter_text>So the value for `2001-01-01` (-867) is now labeled `2000-02-01`. We can collect many of these with a list comprehension and a `concat`.<jupyter_code>lagged = pd.concat([y_prime.shift(i) for i in range(9)], axis=1,
                   keys=['y', 'L1', 'L2', 'L3', 'L4', 'L5', 'L6', 'L7', 'L8'])
lagged
mod_lagged = smf.ols('y ~ L1 + L2 + L3 + L4 + L5 + L6 + L7 + L8', lagged)
res_lagged = mod_lagged.fit()

res_lagged.summary()
ax = res_lagged.fittedvalues.plot(label="predicted", figsize=(12, 4), legend=True)
y_prime.plot(label="actual", legend=True);<jupyter_output><empty_output><jupyter_text>In practice, you won't be doing the `shift`ing and `diff`ing yourself.It's more convenient to let statsmodels do that for us.Then we don't have to worry about un-differencing the fitted / predicted results to interpret them correctly.Also, the solvers we'll see next are a bit more sophisticated than a linear regression.  AutoRegressive ModelPredict $y_{t+1}$, given $y_0, y_1, \ldots y_t$ Let's fit an autoregressive (AR) model. Autoregressive part just means using past values of $y$ to predict the future (like we did above).We'll use statsmodel's `SARIMAX` model. The AR part of SARIMAX is for autoregressive.It also handles seasonality (**S**), differencing (**I** for integrated), moving average (**MA**), and exogenous regressors (**X**).We'll stick to a simple AR(8) model (use the last 8 periods) with a single period of differencing.<jupyter_code>mod = smt.SARIMAX(y, order=(8, 1, 0))  # AR(8), first difference, no MA
res = mod.fit()<jupyter_output><empty_output><jupyter_text>As usual with statsmodels, we get a nice summary with the fitted coefficeints and some test statistics (which we'll ignore)<jupyter_code>res.summary()<jupyter_output><empty_output><jupyter_text>The results instance has all the usual attributes and methods, like `fittedvalues`.<jupyter_code>ax = res.fittedvalues.iloc[1:].plot(label="Fitted", legend=True, figsize=(12, 4))
y.plot(ax=ax, label="Actual", legend=True);<jupyter_output><empty_output><jupyter_text>ForecastingThe real value of timeseries analysis is to predict the future.We can use the `.get_prediction` method to get the predicted values, along with a confidence interval. First, we'll look at one-period-ahead forecasts.Basically, this simulates looking at our data the last day of the month, and making the forecast for the next month.Keep in mind though that we fit our parameters on the entire dataset. The isn't an out-of-sample prediction.<jupyter_code>pred = res.get_prediction(start='2001-03-01')
pred_ci = pred.conf_int()
ax = y.plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='Forecast', alpha=.7)
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)
plt.legend()
sns.despine()<jupyter_output><empty_output><jupyter_text>Alternatively, we can make dynamic forecasts as of some month (January 2013 in the example below). That means the forecast from that point forward only use information available as of January 2013 (though again, we fit the model on the entire dataset). The predictions are generated in a similar way: a bunch of one-step forecasts. Only instead of plugging in the actual values beyond January 2013, we plug in the forecast values.<jupyter_code>pred_dy = res.get_prediction(start='2002-03-01', dynamic='2013-01-01')
pred_dy_ci = pred_dy.conf_int()
ax = y.plot(label='observed')
pred_dy.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_dy_ci.index,
                pred_dy_ci.iloc[:, 0],
                pred_dy_ci.iloc[:, 1], color='k', alpha=.25)
ylim = ax.get_ylim()
ax.fill_betweenx(ylim, pd.Timestamp('2013-01-01'), y.index[-1],
                 alpha=.1, zorder=-1)
ax.set_ylim(ylim)
ax.annotate('Dynamic $\\longrightarrow$',
            (pd.Timestamp('2013-02-01'), 16000))

plt.legend()
sns.despine()
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10927914110429447, 'fraction_numerical': 0.01572085889570552, 'mean_word_length': 5.842077649527807, 'pattern_counts': {'":': 0, '<': 129, '<?xml version=': 0, '>': 129, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '62508', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4348, 'n_tokens_neox': 4082, 'n_words': 1321, 'path': 'notebooks/05-Timeseries.ipynb', 'repo_name': 'eseiver/Pandas-Tutorial-SciPyConf-2018', 'size': '13040'}