%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
import seaborn as sns

# read data into a pandas dataframe
df = pd.read_csv('NH_seaice_extent_final_v2.csv')

# light munging to deal with file format and dump unnecessary columns
df.columns = [col.strip() for col in df.columns.tolist()]
df = df.drop(0) # this is a descriptor row, not actual data
df = df.drop(['Missing','Source Data'], axis=1)
df.Extent = pd.to_numeric(df.Extent)
df.Year = pd.to_numeric(df.Year)
df.Month = pd.to_numeric(df.Month)
df.head()
# Year, Month, Day, Extent (in 10^6 square km)

monthly = df.groupby(['Year', 'Month'])['Extent'].mean()
monthly = monthly.reset_index()
#monthly = monthly.reset_index(level=0)
monthly['Day'] = 1
monthly.index = pd.to_datetime(monthly.loc[:, ['Year', 'Month', 'Day']])

monthly.head()

"""
Explanation: An Intro to Time Series Analysis in Python
The goal of this notebook is introduce some key concepts for getting started with time series analysis. After working through this short notebook, you should have a sense of a where to start with analyzing time series data and how you can begin generating insights.
This tutorial leans heavily on content found in here and here, both of which can be used as a jumping off point for reading.  While useful, these are much longer tutorials. Our goal here is to provide you with something streamlined to get you started quickly!
Note: Please send along any errata or comments you may have.
Outline


Loading example data
The Concept of Stationarity
Assessing Stationarity
Looking at the data
Augmented Dickey-Fuller Test
<br>
<br>


Identifying Trend and Seasonality
Removing Trend to Make a Time Series Stationary
Differencing to Detrend     
Decomposition - Modeling Trend and Season Separately   
<br>


Relating Two Separate Time Series
Granger Causality
<br>
<br>


Take Aways and Next Steps

Dependencies

This notebook uses data from the files "NH_seaice_extent_final_v2.csv" and "Antarctic_Temperatures.csv". In addition, this notebook requires the following libraries:
- numpy
- pandas
- matplotlib
- seaborn
- statsmodels
<a id='data'></a>
Loading Example Data
The data we will work with in this tutorial is real measurements of Antarctic Ice extent. As you can imagine, there is probably a strong seasonal aspect to ice melting and re-freezing, as well as possible long term climate change trends. The original data were on a daily timescale, but we will average within months to simplify things a bit.
You should be able to just run these lines to get your data prepped for the tutorial
End of explanation
"""


# plot ice extent over time
sns.set_context('talk')
monthly.plot(x=monthly.index, y='Extent')
sns.despine()

"""
Explanation: <a id='stationarity-concept'></a>
The Concept of Stationary Data
Stationary Timeseries data have constant mean and variance over time, as well as constant autocorrelation over time. If you look at the seasonality, the frequency and amplitude should remain somewhat constant over time. Additionally, there shouldn't be a trend for the mean to increase or decrease as time goes on.


Of course, real data doesn't behave like this, it almost always is non-stationary. In this example, we'll look at data describing the extent of antarctic ice measured monthly. You can imagine that this varies both within a year (a season component), as well as across years (a trend).
A good discussion of the importance of stationary data to time series analysis can be found here. A key take away on the importance of stationarity is that it is a necessary condition for many of the assumptions that underlie statistical tests.  In short, many of the assumptions we can apply to independent random variables also apply to stationary time series data.
For almost all time series analyses we are going to want to apply (at least some) statistical tests before we dive into modeling.  Therefore, the first thing to do when working with time series data is to assess whether it is stationary.
<a id='stationarity-assessing'></a>
Assessing if Data Are Stationary
<a id='eyeballs'></a>
Looking at the Data
You can often simply plot data and eyeball whether or not they are stationary. In the plot below, there is a small decreasing trend, but it is hard to tell whether or not the time series is actually stationary.
End of explanation
"""


# calculate rolling mean: win = size of window to calculate over
# the window should usually be at least the size of your seasonal period (12 months here)
win = 12
rollMean = monthly.Extent.rolling(window=win, center=False).mean()
rollStd = monthly.Extent.rolling(window=win, center=False).std()

# plot the rolling stats
sns.set_context('talk')
fig, (ax1, ax2) = plt.subplots(2)
ax1.plot(monthly.Extent, label='Raw Data')
ax1.plot(rollMean, label='Rolling Mean')
ax1.set_ylabel('Mean Extent')
ax1.legend()
ax2.plot(rollStd, label='Rolling Std')
ax2.set_xlabel('Year')
ax2.set_ylabel('Extent Std')
sns.despine()

"""
Explanation: One useful tool can be to plot the rolling mean, help you to visualize any trends over time. Additionally, it is smart to plot the standard deviation as a line graph to see if it deviates from a straight line.
End of explanation
"""


from statsmodels.tsa.stattools import adfuller

dftest = adfuller(monthly.Extent, autolag='AIC')

# the statsmodels output is ugly and unlabeled, lets use a function to format it nicely
def adf_output(dftest):
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','# of Lags Used','Number of Observations Used'])

    # the last element is a dictionary to unpack
    for key,value in dftest[4].items(): 
        dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)
    
adf_output(dftest)

"""
Explanation: Based on these plots, it looks like neither the mean nor the standard deviation are stable over time. But it is still only a small variation, so it is still hard to judge whether or not this time series is stationary. Instead, we can use statistical tests to say if the data are stationary.
<a id='adf'></a>
Augmented Dickey-Fuller Test
Time Series models are described in terms of autoregressive (AR) lags. So an AR(1) model states that the value at Y<sub>t</sub> depends on the value at Y<sub>t-1</sub>. So we can represent the time series as a model:
Y<sub>t</sub> = Alpha + Rho*Y<sub>t-1</sub> + Error<sub>t</sub>
It turns out, a time series is perfectly stationary when Rho = 0. This makes intuitive sense: The time series will not depend on its previous value when Rho is 0, so it depends entirely on error and is equivalent to a random walk. A time series is perfectly non-stationary if Rho = 1. Alpha is a coeffecient that we want to estimate such that the Error term is minimized, but let's worry about that later.
The Dickey-Fuller Test asks if Rho in an AR(1) model is equal to 0 (H<sub>0</sub>: Rho=1, H<sub>1</sub>: Rho<1)
The augmented Dickey-Fuller Test just expands this to work for AR(n) models
A p-value < 0.05 indicates that your time series IS STATIONARY.
End of explanation
"""


# calculate rolling mean
win = 12
rollMean = monthly.Extent.rolling(window=win, center=False).mean()

# plot
sns.set_context('talk')
plt.plot(monthly.Extent, label='Raw Data')
plt.plot(rollMean, color='red', label='Rolling Mean')
plt.xlabel('Year')
plt.ylabel('Extent')
plt.legend()
plt.title("Raw Data and the Rolling Mean")
sns.despine()

"""
Explanation: In this output we can see that the p-value is much larger that 0.05, confirming that the data are not stationary!
<a id='identifying-trends'></a>
Identifying Trend and Seasonality
Common sources of non-stationarity are Trend and Seasonality. Trend is when your mean drifts over time, while seasonality is variations in specific time frames (like ice melting in summer but accumulating in winter).
In addition to removing these components for forecasting purposes, most data scientists will want to identify these components and look at them historically to see if there are possible reasons for changes in them over time. Perhaps the trend in ice melting reverses around the time electric cars are introduced, for example.
<a id='trend-removal'></a>
Removing Trend to Make a Time Series Stationary
The simplest way to estimate and remove trend is to calculate a rolling mean and subtract it out:
End of explanation
"""


# subtract mean
demeaned = monthly.Extent - rollMean
demeaned.dropna(inplace=True)

# did this do anything for our stationarity?
demeaned_adf = adfuller(demeaned, autolag='AIC')
adf_output(demeaned_adf)

"""
Explanation: Notice that in the above plot, a few data points are missing from the front of the rolling mean time series. This is because of the window for calculating the moving average. We can remove those NAs and then detrend the time series. This loss of data is a good reason to keep the window as small as possible.
End of explanation
"""


# rolling mean and std
win = 12
roll_demean = demeaned.rolling(window=win, center=False).mean()
roll_demean_std = demeaned.rolling(window=win, center=False).std()

# plot
plt.plot(roll_demean, color='red', label='mean')
plt.plot(roll_demean_std, color='black', label='std')
plt.legend(loc=2)
sns.despine()

"""
Explanation: With the output above, you can see that this simple transformation made the data stationary. <br> Below, you can see that the moving average is now more stable.
End of explanation
"""


# plot the de-meaned time series
plt.plot(demeaned)
sns.despine()

"""
Explanation: This is also apparent in the lack of decreasing trend in our demeaned time series
End of explanation
"""


# construct difference curve and plot
diffed = monthly.Extent - monthly.Extent.shift()
diffed.dropna(inplace=True)
plt.plot(diffed)
sns.despine()

diffed_adf = adfuller(diffed, autolag='AIC')
adf_output(diffed_adf)

"""
Explanation: <a id='differencing'></a>
Differencing to Detrend
Another common method for removing trend is differencing. In this technique, you just subtract from your time series the time series shifted by t-1. Much like demeaning, this operation leaves you with just the change across time, which  often leads to a stationary series. You can then just model that change and superimpose it back on the trend.
End of explanation
"""


from statsmodels.tsa.seasonal import seasonal_decompose

decomp = seasonal_decompose(monthly.Extent, model="additive")

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4)
ax1.plot(monthly.Extent)
ax2.plot(decomp.trend)
ax3.plot(decomp.seasonal)
ax4.plot(decomp.resid)
sns.despine()

"""
Explanation: Differencing can also be useful when done multiple times, depending on the time series.
<br>
<a id='decomposition'></a>
Decomposition - Modeling Trend and Season Separately
End of explanation
"""


resid = decomp.resid.dropna()
decomp_adf = adfuller(resid, autolag='AIC')
adf_output(decomp_adf)

"""
Explanation: It is worth noting that statsmodels themselves admit this is a naive decomposition. The R statistical language has substantially more support for time series analysis than does Python, so it may be worth exploring STL Decomposition in R if this decomposition is not getting you very far.
Regardless of how you do your decomposition, the residual component is what you would use to do your modeling.
A major benefit of decomposition is that now you can visualize the trend and seasonal components separately, allowing you to make inferences from those sources of data.
End of explanation
"""


# add monthly temperature data to the mix
tmp = pd.read_csv("Antarctic_Temperatures.csv")
av_temp = pd.melt(tmp, id_vars=['Year'], value_vars=['1','2','3','4','5','6','7','8','9','10','11','12'])
av_temp.columns = ['Year', 'Month', 'Temperature']
av_temp['Day'] = 1
av_temp.index = pd.to_datetime(av_temp.ix[:,['Year', 'Month', 'Day']])
av_temp.head()

monthly['Temperature'] = av_temp.Temperature
monthly.head()

plt.plot((monthly.Extent - np.mean(monthly.Extent)) / np.std(monthly.Extent))
plt.plot((monthly.Temperature - np.mean(monthly.Temperature)) / np.std(monthly.Temperature))
sns.despine()

"""
Explanation: <a id='two-series'></a>
Relating Two Separate Time Series
You might have some other time series and want to see if it is related to changes in your first time series. For example, we might expect that monthly temperatures predict ice extent in Antarctica
End of explanation
"""


temp_adf = adfuller(monthly.Temperature, autolag='AIC')
adf_output(temp_adf)

"""
Explanation: While hard to see in this graph, it looks like peaks in temperature might precede drops in ice coverage. Lets statistically test this relationship
<a id='granger'></a>
Granger Causality
Granger Causality is a statistical method for testing the relationship between two time series. This can be a great way to supplement an analysis and move beyond simply plotting the decomposed time series.
At it's heart, this method asks whether values from a second time series provide more information about the future values of your first time series than just considering past values of the first time series:
"X is said to Granger-cause Y if Y can be better predicted using the histories of both X and Y than it can by using the history of Y alone." A good explanation of Granger Causality can be found here and here (the first is much simpler, the second is for a real deep dive!).
An important caveat of Granger Causality is that it is very possible to find that X significantly Granger-causes Y, and simultaneously see that Y Granger-causes X. With well selected variables, you can curb this issue to some degree. For example, while maybe possible, it is much more likely that changes in temperature cause ice melting than ice melting causing changes in temperature. It may also be easier to get unidirectional causality in trend components rather than seasonal components.
 First, both time series should be stationary. So we'll do an ADF test on the temperature data and then use the differenced data for monthly ice extent
End of explanation
"""


from statsmodels.tsa.stattools import grangercausalitytests

# create a 2d dataframe to test if ice extent is caused by temperature
granger = pd.DataFrame({'diffed_extent':diffed, 'temperature':av_temp.Temperature})
granger = granger.dropna()

granger_test = grangercausalitytests(granger, maxlag=12)

"""
Explanation: The Statsmodels Granger Causality package is a bit hard to parse: 
- You should provide it a 2 dimensional structure, where the first column is expected to be caused by the second column
- the maxlag argument lets you run the test with lagged data from your second column
- statsmodels will give you results for a test, incrementing lags, all the way up to your max lags
End of explanation
"""


