<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Time series basics
Trends
Seasonality
Cyclical
Introduction to statsmodels
Step1:  Using the Hodrick-Prescott Filter for trend analysis
Step2:  ETS Theory (Error-Trend-Seasonality)
Exponential Smoothing
Trend Methods Models
ETS Decomposition
EWMA Theory
(Exponentially Weighted Moving Averages)
Step3:  Weakness of SMA
* Smaller windows will lead to more noise, rather than signal
* It will always lag by the size of the window
* It will never reach to full peak or valley of the data due to the averaging.
* Does not really inform you about possible future behaviour, all it really does is describe trends in your data.
* Extreme historical values can skew your SMA significantly
Creating EWMA
Step4:  Full reading on mathematics of EWMA
* http
Step5:  ARIMA models
Auto regressive integrated moving averages
[https
Step6:  Step 1 - Visualize the data
Step7:  Conclusion
Step8:  ARIMA models continued... 2
### Step 2 - Make the time series data stationary (if non-stationary)
We can use the Augmented Dickey-Fuller unit root test.
In statistics and econometrics, an augmented Dickey–Fuller test (ADF) tests the null hypothesis that a unit root is present in a time series sample. The alternative hypothesis is different depending on which version of the test is used, but is usually stationarity or trend-stationarity.
Basically, we are trying to whether to accept the Null Hypothesis H0 (that the time series has a unit root, indicating it is non-stationary) or reject H0 and go with the Alternative Hypothesis (that the time series has no unit root and is stationary).
We end up deciding this based on the p-value return.
A small p-value (typically ≤ 0.05) indicates strong evidence against the null hypothesis, so you reject the null hypothesis.
A large p-value (> 0.05) indicates weak evidence against the null hypothesis, so you fail to reject the null hypothesis.
Let's run the Augmented Dickey-Fuller test on our data
Step9:  Thus, the ADF test confirms our assumption from visual analysis that definately the data is non-stationary and has a seasionality and trend factor to it.
Step10:  Since, p-value('Original') - p-value('First Difference') < p-value('First Difference') - p-value('Second Difference'), it is the first difference that did most of the elimination of the trend.
Step11:  Thus, we conclude that seasonal difference does not make the data stationary here, in fact we can observe visually that as we go further in time the variance began to increase.
Step12:  ARIMA models continued... 3
### Step 3 - Plot the Correlation and Autocorrelation Charts
Step13:  Plotting the final 'Autocorrelation' and 'Partial autocorrelation'
Step14:  ARIMA models continued... 4
### Step 4 - Construct the ARIMA model
Step15:  Choosing the p, d, q values of the order and seasonal_order tuple is reading task
More information here...
[https
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import statsmodels.api as sm
# Importing built-in datasets in statsmodels
df = sm.datasets.macrodata.load_pandas().data
df.head()
print(sm.datasets.macrodata.NOTE)
df.head()
df.tail()
# statsmodels.timeseriesanalysis.datetools
index = pd.Index(sm.tsa.datetools.dates_from_range('1959Q1', '2009Q3'))
index
df.index = index
df.head()
df['realgdp'].plot()
Explanation: Time series basics
Trends
Seasonality
Cyclical
Introduction to statsmodels
End of explanation
result = sm.tsa.filters.hpfilter(df['realgdp'])
result
type(result)
type(result[0])
type(result[1])
gdp_cycle, gdp_trend = result
df['trend'] = gdp_trend
df[['realgdp', 'trend']].plot()
# zooming in
df[['realgdp', 'trend']]['2000-03-31':].plot()
Explanation: Using the Hodrick-Prescott Filter for trend analysis
End of explanation
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
airline = pd.read_csv('airline_passengers.csv', index_col = 'Month')
airline.head()
# this is a normal index
airline.index
# Get rid of all the missing values in this dataset
airline.dropna(inplace=True)
airline.index = pd.to_datetime(airline.index)
airline.head()
# now its a DatetimeIndex
airline.index
# Recap of making the SMA
airline['6-month-SMA'] = airline['Thousands of Passengers'].rolling(window=6).mean()
airline['12-month-SMA'] = airline['Thousands of Passengers'].rolling(window=12).mean()
airline.plot(figsize=(10,8))
Explanation: ETS Theory (Error-Trend-Seasonality)
Exponential Smoothing
Trend Methods Models
ETS Decomposition
EWMA Theory
(Exponentially Weighted Moving Averages)
End of explanation
airline['EWMA-12'] = airline['Thousands of Passengers'].ewm(span=12).mean()
airline[['Thousands of Passengers', 'EWMA-12']].plot(figsize=(10,8))
Explanation: Weakness of SMA
* Smaller windows will lead to more noise, rather than signal
* It will always lag by the size of the window
* It will never reach to full peak or valley of the data due to the averaging.
* Does not really inform you about possible future behaviour, all it really does is describe trends in your data.
* Extreme historical values can skew your SMA significantly
Creating EWMA
End of explanation
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
airline = pd.read_csv('airline_passengers.csv', index_col = 'Month')
airline.head()
airline.plot()
airline.dropna(inplace = True)
airline.index = pd.to_datetime(airline.index)
airline.head()
from statsmodels.tsa.seasonal import seasonal_decompose
# additive/ multiplicative models available
# suspected linear trend = use additive
# suspected non-linear trend = multiplicative model
result = seasonal_decompose(airline['Thousands of Passengers'], model='multiplicative')
result.seasonal
result.seasonal.plot()
result.trend.plot()
result.resid.plot()
fig = result.plot()
Explanation: Full reading on mathematics of EWMA
* http://pandas.pydata.org/pandas-docs/stable/computation.html#exponentially-weighted-windows
ETS continued...
End of explanation
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline
df = pd.read_csv('monthly-milk-production-pounds-p.csv')
df.head()
df.columns = ['Month', 'Milk in Pounds per Cow']
df.head()
df.tail()
df.drop(168, axis=0, inplace=True)
df.tail()
df['Month'] = pd.to_datetime(df['Month'])
df.head()
df.set_index('Month', inplace=True)
df.head()
df.index
df.describe()
df.describe().transpose()
Explanation: ARIMA models
Auto regressive integrated moving averages
[https://people.duke.edu/~rnau/411arim3.htm]
Autoregressive integrated moving average (ARIMA) model is a generalization of an autoregressive moving average (ARMA) model.
ARIMA model types
Non-seasonal ARIMA (for non-seasonal data)
Seasonal ARIMA (for seasonal data)
ARIMA models are applied in some cases where data show evidence of non-stationarity, where an initial differencing step (corresponding to the 'integrated' part of the model) can be applied one or more times to eliminate the non-stationarity.
Non-seasonal ARIMA models are generally denoted as ARIMA(p, d, q) where parameters p, d and q are non-negative integers.
AR(p): Autoregression component
A regression model that utilizes the dependent relationship between a current observation and observations over a previous period.
I(d): Integrated
Differencing of observations (subtracting an observation from an observation at the previous time step) in order to make the time series stationary.
MA(q): Moving Average
A model that uses the dependency between an observation and a residual error from a moving average model applied to lagged observations
Stationary vs Non-Stationary Data
A stationary series has a constant mean and variance over time
A stationary dataset will allow our model to predict that the mean and variance will be the same in future periods.
<img src='stationaryvsnonstationary.png' />
Note above for stationary data (mean and variance both are constant over time)
Another aspect to look for is covariance not be a function of time in stationary data
If you've determined your data is not stationary (either visually or mathematically), you will then need to transform it to be stationary in order to evaluate it and what type of ARIMA terms you will use.
One simple way to do this is through "differencing"
Original Data
<table>
<tr><td>Time1</td><td>10</td></tr>
<tr><td>Time2</td><td>12</td></tr>
<tr><td>Time3</td><td>8</td></tr>
<tr><td>Time4</td><td>14</td></tr>
<tr><td>Time5</td><td>7</td></tr>
</table>
First Difference
<table>
<tr><td>Time1</td><td>NA</td></tr>
<tr><td>Time2</td><td>2</td></tr>
<tr><td>Time3</td><td>-4</td></tr>
<tr><td>Time4</td><td>6</td></tr>
<tr><td>Time5</td><td>-7</td></tr>
</table>
Second Difference
<table>
<tr><td>Time1</td><td>NA</td></tr>
<tr><td>Time2</td><td>NA</td></tr>
<tr><td>Time3</td><td>-6</td></tr>
<tr><td>Time4</td><td>10</td></tr>
<tr><td>Time5</td><td>-13</td></tr>
</table>
For seasonal data, you can also difference by season. If you had monthly data with yearly seasonality, you could difference by a time unit of 12, instead of just 1.
Another common techinique with seasonal ARIMA models is to combine both methods, taking the seasonal difference of the first difference.
ARIMA models continued... 1
The general process for ARIMA models is the following:
* Visualize the Time Series Data
* Make the time series data stationary
* Plot the Correlation and AutoCorrelation Charts
* Construct the ARIMA Model
* Use the model to make predictions
Let's go through these steps!
[https://people.duke.edu/~rnau/arimrule.htm]
End of explanation
df.plot();
time_series = df['Milk in Pounds per Cow']
type(time_series)
time_series.rolling(12).mean().plot(label='12 SMA')
time_series.rolling(12).std().plot(label='12 STD')
time_series.plot()
plt.legend();
Explanation: Step 1 - Visualize the data
End of explanation
from statsmodels.tsa.seasonal import seasonal_decompose
decomp = seasonal_decompose(time_series)
fig = decomp.plot()
fig.set_size_inches(15,8)
Explanation: Conclusion: The scale of STD (standard deviation) is always pretty much smaller than the actual scale. If 12 STD does not show crazy behaviour is comparitively flat then its 'workable'
End of explanation
from statsmodels.tsa.stattools import adfuller
result = adfuller(df['Milk in Pounds per Cow'])
result
def adf_check(time_series):
    result = adfuller(time_series)
    print('Augumented Dicky-Fuller Test')
    labels = ['ADF Test Statistic', 'p-value', '# of lags', 'Num of Observations used']
    for value, label in zip(result, labels):
        print(label + ': ' + str(value))
        
    if result[1] < 0.05:
        print('Strong evidence against null hypothesis')
        print('Rejecting null hypothesis')
        print('Data has no unit root! and is stationary')
    else:
        print('Weak evidence against null hypothesis')
        print('Fail to reject null hypothesis')
        print('Data has a unit root, it is non-stationary')
adf_check(df['Milk in Pounds per Cow'])
Explanation: ARIMA models continued... 2
### Step 2 - Make the time series data stationary (if non-stationary)
We can use the Augmented Dickey-Fuller unit root test.
In statistics and econometrics, an augmented Dickey–Fuller test (ADF) tests the null hypothesis that a unit root is present in a time series sample. The alternative hypothesis is different depending on which version of the test is used, but is usually stationarity or trend-stationarity.
Basically, we are trying to whether to accept the Null Hypothesis H0 (that the time series has a unit root, indicating it is non-stationary) or reject H0 and go with the Alternative Hypothesis (that the time series has no unit root and is stationary).
We end up deciding this based on the p-value return.
A small p-value (typically ≤ 0.05) indicates strong evidence against the null hypothesis, so you reject the null hypothesis.
A large p-value (> 0.05) indicates weak evidence against the null hypothesis, so you fail to reject the null hypothesis.
Let's run the Augmented Dickey-Fuller test on our data:
End of explanation
# Now making the data stationary
df['First Difference'] = df['Milk in Pounds per Cow'] - df['Milk in Pounds per Cow'].shift(1)
df['First Difference'].plot()
# adf_check(df['First Difference']) - THIS RESULTS IN LinAlgError: SVD did not converge ERROR
# Note: we need to drop the first NA value before plotting this
adf_check(df['First Difference'].dropna())
df['Second Difference'] = df['First Difference'] - df['First Difference'].shift(1)
df['Second Difference'].plot();
adf_check(df['Second Difference'].dropna())
Explanation: Thus, the ADF test confirms our assumption from visual analysis that definately the data is non-stationary and has a seasionality and trend factor to it.
End of explanation
# Let's plot seasonal difference
df['Seasonal Difference'] = df['Milk in Pounds per Cow'] - df['Milk in Pounds per Cow'].shift(12)
df['Seasonal Difference'].plot();
adf_check(df['Seasonal Difference'].dropna())
Explanation: Since, p-value('Original') - p-value('First Difference') < p-value('First Difference') - p-value('Second Difference'), it is the first difference that did most of the elimination of the trend.
End of explanation
# Plotting 'Seasonal first difference'
df['Seasonal First Difference'] = df['First Difference'] - df['First Difference'].shift(12)
df['Seasonal First Difference'].plot();
adf_check(df['Seasonal First Difference'].dropna())
Explanation: Thus, we conclude that seasonal difference does not make the data stationary here, in fact we can observe visually that as we go further in time the variance began to increase.
End of explanation
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# Plotting the gradual decline autocorrelation
fig_first = plot_acf(df['First Difference'].dropna())
fig_first = plot_acf(df['First Difference'].dropna(), use_vlines=False)
fig_seasonal_first = plot_acf(df['Seasonal First Difference'].dropna(), use_vlines=False)
fig_seasonal_first_pacf = plot_pacf(df['Seasonal First Difference'].dropna(), use_vlines=False)
Explanation: ARIMA models continued... 3
### Step 3 - Plot the Correlation and Autocorrelation Charts
End of explanation
plot_acf(df['Seasonal First Difference'].dropna());
plot_pacf(df['Seasonal First Difference'].dropna());
Explanation: Plotting the final 'Autocorrelation' and 'Partial autocorrelation'
End of explanation
# ARIMA model for non-sesonal data
from statsmodels.tsa.arima_model import ARIMA
# help(ARIMA)
# ARIMA model from seasonal data
# from statsmodels.tsa.statespace import sarimax
Explanation: ARIMA models continued... 4
### Step 4 - Construct the ARIMA model
End of explanation
model = sm.tsa.statespace.SARIMAX(df['Milk in Pounds per Cow'], order=(0,1,0), seasonal_order=(1,1,1,12))
results = model.fit()
print(results.summary())
# residual errors of prediction on the original training data
results.resid
# plot of residual errors of prediction on the original training data
results.resid.plot();
# KDE plot of residual errors of prediction on the original training data
results.resid.plot(kind='kde');
# Creating a column forecast to house the forecasted values for existing values
df['forecast'] = results.predict(start=150, end=168)
df[['Milk in Pounds per Cow', 'forecast']].plot(figsize=(12,8));
# Forecasting for future data
df.tail()
from pandas.tseries.offsets import DateOffset
future_dates = [df.index[-1] + DateOffset(months=x) for x in range(0,24)]
future_dates
future_df = pd.DataFrame(index=future_dates, columns=df.columns)
future_df.head()
final_df = pd.concat([df, future_df])
final_df.head()
final_df.tail()
final_df['forecast'] = results.predict(start=168, end=192)
final_df.tail()
final_df[['Milk in Pounds per Cow', 'forecast']].plot()
Explanation: Choosing the p, d, q values of the order and seasonal_order tuple is reading task
More information here...
[https://stackoverflow.com/questions/22770352/auto-arima-equivalent-for-python]
[https://stats.stackexchange.com/questions/44992/what-are-the-values-p-d-q-in-arima]
[https://people.duke.edu/~rnau/arimrule.htm]
End of explanation
<END_TASK>