In [1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from StringIO import StringIO
import requests
print 'Numpy:', np.__version__
print 'Pandas:', pd.__version__
import statsmodels
print 'Statsmodels:', statsmodels.__version__
np.set_printoptions(precision=2)
In [2]:
act = requests.get('https://docs.google.com/spreadsheet/ccc?key=0Ak_wF7ZGeMmHdFZtQjI1a1hhUWR2UExCa2E4MFhiWWc&output=csv&gid=1')
data = pd.read_csv(StringIO(act.content),index_col=0,parse_dates=['date'], thousands=',').sort_index() #converts to numbers
In [3]:
# In your example, you weren't able to predict because your dataframe was not sorted correctly;
# If you look at the model summary table, you'll see that the sample runs from
# 08-01-2013 to 05-01-2000
data['Unit Sales'].plot()
Out[3]:
In [4]:
# Augmented Dicky-Fuller test for unit roots
# Null Hypothesis: Unit Root
# Alternative Hypothesis: AR(p), where here p is selected as the lag <= 12 that minimizes AIC
print 'Augmented Dicky-Fuller Test (p <= 12)'
print sm.tsa.adfuller(data['Unit Sales'], 12)[1], '=> we cannot reject unit root in data'
print sm.tsa.adfuller(data['Unit Sales'].diff()[1:], 11)[1], '=> we can reject unit root in differenced data at 1%'
In [5]:
# So it looks like an appropriate model might be ARIMA(12,1,0)
p=12
d=1
q=0
mod = sm.tsa.ARIMA(data['Unit Sales'], order=[p,d,q], freq='M') # mod is an object of class ARIMA
res = mod.fit(); # res is an object of class ARIMAResults
In [6]:
print res.summary() # sometimes the Display() command isn't as nice as print
# The roots table gives you information about the stationarity of the estimated ARIMA
# model (i.e. did you choose a high enough `d`). All the roots are greater than 1 in
# modulus, so our model is stationary. We could guess this would happen based on the
# ADF test of the differenced data, which told us that we had one unit root.
In [7]:
predict = res.predict('2012-1-01','2015-01-01', dynamic=True)
fig = plt.figure(figsize=(10,3))
ax = fig.add_subplot(111)
#ax.plot(data.index, data['Unit Sales'])
ax.plot(data.index, data['Unit Sales']-data['Unit Sales'].mean())
ax.plot(predict.index, predict, 'r-');
In [7]: