Autoregressive Moving Average (ARMA): Sunspots data

This notebook replicates the existing ARMA notebook using the statsmodels.tsa.statespace.SARIMAX class rather than the statsmodels.tsa.ARMA class.

In [1]:
%matplotlib inline
In [2]:
from __future__ import print_function
import numpy as np
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt

import statsmodels.api as sm
/build/statsmodels-GE7Zhw/statsmodels-0.8.0/.pybuild/cpython3_3.6_statsmodels/build/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools
In [3]:
from statsmodels.graphics.api import qqplot

Sunpots Data

In [4]:
print(sm.datasets.sunspots.NOTE)
::

    Number of Observations - 309 (Annual 1700 - 2008)
    Number of Variables - 1
    Variable name definitions::

        SUNACTIVITY - Number of sunspots for each year

    The data file contains a 'YEAR' variable that is not returned by load.

In [5]:
dta = sm.datasets.sunspots.load_pandas().data
In [6]:
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))
del dta["YEAR"]
In [7]:
dta.plot(figsize=(12,4));
In [8]:
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(dta.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(dta, lags=40, ax=ax2)
In [9]:
arma_mod20 = sm.tsa.statespace.SARIMAX(dta, order=(2,0,0), trend='c').fit()
print(arma_mod20.params)
intercept     14.793947
ar.L1          1.390659
ar.L2         -0.688568
sigma2       274.761104
dtype: float64
In [10]:
arma_mod30 = sm.tsa.statespace.SARIMAX(dta, order=(3,0,0), trend='c').fit()
In [11]:
print(arma_mod20.aic, arma_mod20.bic, arma_mod20.hqic)
2622.6363381416013 2637.569703249192 2628.6067259868473
In [12]:
print(arma_mod30.params)
intercept     16.762205
ar.L1          1.300810
ar.L2         -0.508122
ar.L3         -0.129612
sigma2       270.102652
dtype: float64
In [13]:
print(arma_mod30.aic, arma_mod30.bic, arma_mod30.hqic)
2619.4036296630375 2638.0703360475263 2626.866614469595
  • Does our model obey the theory?
In [14]:
sm.stats.durbin_watson(arma_mod30.resid)
Out[14]:
1.9564844846369152
In [15]:
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
ax = plt.plot(arma_mod30.resid)
In [16]:
resid = arma_mod30.resid
In [17]:
stats.normaltest(resid)
Out[17]:
NormaltestResult(statistic=49.847005587191916, pvalue=1.499202393979229e-11)
In [18]:
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
In [19]:
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
In [20]:
r,q,p = sm.tsa.acf(resid, qstat=True)
data = np.c_[range(1,41), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
print(table.set_index('lag'))
            AC          Q      Prob(>Q)
lag                                    
1.0   0.009176   0.026273  8.712350e-01
2.0   0.041820   0.573727  7.506142e-01
3.0  -0.001342   0.574292  9.022915e-01
4.0   0.136064   6.407489  1.707135e-01
5.0   0.092433   9.108335  1.048203e-01
6.0   0.091919  11.788019  6.686839e-02
7.0   0.068735  13.291376  6.531938e-02
8.0  -0.015021  13.363413  9.994244e-02
9.0   0.187599  24.636917  3.400196e-03
10.0  0.213724  39.317882  2.233181e-05
11.0  0.201092  52.358271  2.347758e-07
12.0  0.117192  56.802110  8.581663e-08
13.0 -0.014051  56.866211  1.895533e-07
14.0  0.015394  56.943404  4.001104e-07
15.0 -0.024986  57.147465  7.747082e-07
16.0  0.080892  59.293627  6.880517e-07
17.0  0.041120  59.850086  1.112485e-06
18.0 -0.052030  60.744065  1.550379e-06
19.0  0.062500  62.038495  1.833801e-06
20.0 -0.010292  62.073719  3.385222e-06
21.0  0.074467  63.924063  3.196543e-06
22.0  0.124962  69.152772  8.984831e-07
23.0  0.093170  72.069533  5.802914e-07
24.0 -0.082149  74.345042  4.715785e-07
25.0  0.015689  74.428333  8.294017e-07
26.0 -0.025049  74.641400  1.367992e-06
27.0 -0.125875  80.040874  3.722921e-07
28.0  0.053215  81.009319  4.717355e-07
29.0 -0.038699  81.523324  6.917766e-07
30.0 -0.016896  81.621649  1.151883e-06
31.0 -0.019286  81.750228  1.869202e-06
32.0  0.105001  85.575149  8.927708e-07
33.0  0.040094  86.134872  1.247384e-06
34.0  0.008834  86.162142  2.047606e-06
35.0  0.014588  86.236785  3.263459e-06
36.0 -0.119334  91.249667  1.084187e-06
37.0 -0.036673  91.724838  1.521456e-06
38.0 -0.046204  92.481862  1.937920e-06
39.0 -0.017775  92.594310  2.989369e-06
40.0 -0.006219  92.608126  4.694971e-06
  • This indicates a lack of fit.
  • In-sample dynamic prediction. How good does our model do?
In [21]:
predict_sunspots = arma_mod30.predict(start='1990', end='2012', dynamic=True)
In [22]:
fig, ax = plt.subplots(figsize=(12, 8))
dta.ix['1950':].plot(ax=ax)
predict_sunspots.plot(ax=ax, style='r');
/usr/lib/python3/dist-packages/ipykernel_launcher.py:2: DeprecationWarning: 
.ix is deprecated. Please use
.loc for label based indexing or
.iloc for positional indexing

See the documentation here:
http://pandas.pydata.org/pandas-docs/stable/indexing.html#ix-indexer-is-deprecated
  
In [23]:
def mean_forecast_err(y, yhat):
    return y.sub(yhat).mean()
In [24]:
mean_forecast_err(dta.SUNACTIVITY, predict_sunspots)
Out[24]:
5.6355503127660205