Autoregressive Moving Average (ARMA): Sunspots data

In [1]:
%matplotlib inline

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
The history saving thread hit an unexpected error (OperationalError('disk I/O error')).History will not be written to the database.
/usr/lib/python3/dist-packages/numpy/core/getlimits.py:214: RuntimeWarning: overflow encountered in nextafter
  if hasattr(umath, 'nextafter')  # Missing on some platforms?
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_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 [2]:
from statsmodels.graphics.api import qqplot

Sunpots Data

In [3]:
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 [4]:
dta = sm.datasets.sunspots.load_pandas().data
In [5]:
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))
del dta["YEAR"]
In [6]:
dta.plot(figsize=(12,8));
In [7]:
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 [8]:
arma_mod20 = sm.tsa.ARMA(dta, (2,0)).fit()
print(arma_mod20.params)
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1341: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  out_full[ind] += zi
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1344: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  out = out_full[ind]
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1350: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  zf = out_full[ind]
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:650: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  elif issubdtype(paramsdtype, complex):
const                49.659457
ar.L1.SUNACTIVITY     1.390656
ar.L2.SUNACTIVITY    -0.688571
dtype: float64
In [9]:
arma_mod30 = sm.tsa.ARMA(dta, (3,0)).fit()
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1341: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  out_full[ind] += zi
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1344: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  out = out_full[ind]
/usr/lib/python3/dist-packages/scipy/signal/signaltools.py:1350: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  zf = out_full[ind]
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:650: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  elif issubdtype(paramsdtype, complex):
In [10]:
print(arma_mod20.aic, arma_mod20.bic, arma_mod20.hqic)
2622.6363380640614 2637.5697031716522 2628.6067259093074
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
In [11]:
print(arma_mod30.params)
const                49.750057
ar.L1.SUNACTIVITY     1.300810
ar.L2.SUNACTIVITY    -0.508093
ar.L3.SUNACTIVITY    -0.129649
dtype: float64
In [12]:
print(arma_mod30.aic, arma_mod30.bic, arma_mod30.hqic)
2619.403628698192 2638.0703350826807 2626.8666135047492
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
  • Does our model obey the theory?
In [13]:
sm.stats.durbin_watson(arma_mod30.resid.values)
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:577: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
Out[13]:
1.9564814686097798
In [14]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax = arma_mod30.resid.plot(ax=ax);
In [15]:
resid = arma_mod30.resid
In [16]:
stats.normaltest(resid)
Out[16]:
NormaltestResult(statistic=49.84496922025755, pvalue=1.5007296344442907e-11)
In [17]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
/usr/lib/python3/dist-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval
In [18]:
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
In [19]:
r,q,p = sm.tsa.acf(resid.values.squeeze(), 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.009178   0.026284  8.712081e-01
2.0   0.041793   0.573033  7.508747e-01
3.0  -0.001335   0.573593  9.024502e-01
4.0   0.136089   6.408898  1.706218e-01
5.0   0.092468   9.111798  1.046871e-01
6.0   0.091948  11.793210  6.674426e-02
7.0   0.068748  13.297163  6.519068e-02
8.0  -0.015020  13.369192  9.976247e-02
9.0   0.187592  24.641857  3.393976e-03
10.0  0.213718  39.321922  2.229539e-05
11.0  0.201082  52.361045  2.345040e-07
12.0  0.117181  56.804084  8.574631e-08
13.0 -0.014055  56.868222  1.893982e-07
14.0  0.015398  56.945460  3.997824e-07
15.0 -0.024967  57.149216  7.741781e-07
16.0  0.080916  59.296667  6.872434e-07
17.0  0.041138  59.853637  1.110986e-06
18.0 -0.052021  60.747328  1.548490e-06
19.0  0.062496  62.041591  1.831711e-06
20.0 -0.010302  62.076879  3.381364e-06
21.0  0.074453  63.926549  3.193705e-06
22.0  0.124955  69.154661  8.978709e-07
23.0  0.093162  72.070923  5.800014e-07
24.0 -0.082152  74.346579  4.713198e-07
25.0  0.015695  74.429935  8.289355e-07
26.0 -0.025037  74.642793  1.367335e-06
27.0 -0.125861  80.041039  3.722706e-07
28.0  0.053225  81.009876  4.716448e-07
29.0 -0.038693  81.523701  6.916879e-07
30.0 -0.016904  81.622120  1.151701e-06
31.0 -0.019296  81.750833  1.868830e-06
32.0  0.104990  85.574960  8.928264e-07
33.0  0.040086  86.134461  1.247551e-06
34.0  0.008829  86.161704  2.047894e-06
35.0  0.014588  86.236342  3.263915e-06
36.0 -0.119329  91.248793  1.084490e-06
37.0 -0.036665  91.723760  1.521973e-06
38.0 -0.046193  92.480409  1.938798e-06
39.0 -0.017768  92.592777  2.990776e-06
40.0 -0.006220  92.606600  4.697133e-06
  • This indicates a lack of fit.
  • In-sample dynamic prediction. How good does our model do?
In [20]:
predict_sunspots = arma_mod30.predict('1990', '2012', dynamic=True)
print(predict_sunspots)
1990-12-31    167.047483
1991-12-31    140.993152
1992-12-31     94.859336
1993-12-31     46.861155
1994-12-31     11.242823
1995-12-31     -4.721109
1996-12-31     -1.166797
1997-12-31     16.185746
1998-12-31     39.021906
1999-12-31     59.449897
2000-12-31     72.170200
2001-12-31     75.376888
2002-12-31     70.436606
2003-12-31     60.731762
2004-12-31     50.201979
2005-12-31     42.076196
2006-12-31     38.114429
2007-12-31     38.454756
2008-12-31     41.963906
2009-12-31     46.869367
2010-12-31     51.423343
2011-12-31     54.399813
2012-12-31     55.321802
Freq: A-DEC, dtype: float64
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:577: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
In [21]:
fig, ax = plt.subplots(figsize=(12, 8))
ax = dta.ix['1950':].plot(ax=ax)
fig = arma_mod30.plot_predict('1990', '2012', dynamic=True, ax=ax, plot_insample=False)
/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
  
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:577: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
In [22]:
def mean_forecast_err(y, yhat):
    return y.sub(yhat).mean()
In [23]:
mean_forecast_err(dta.SUNACTIVITY, predict_sunspots)
Out[23]:
5.636825493451873

Exercise: Can you obtain a better fit for the Sunspots model? (Hint: sm.tsa.AR has a method select_order)

Simulated ARMA(4,1): Model Identification is Difficult

In [24]:
from statsmodels.tsa.arima_process import arma_generate_sample, ArmaProcess
In [25]:
np.random.seed(1234)
# include zero-th lag
arparams = np.array([1, .75, -.65, -.55, .9])
maparams = np.array([1, .65])

Let's make sure this model is estimable.

In [26]:
arma_t = ArmaProcess(arparams, maparams)
In [27]:
arma_t.isinvertible
Out[27]:
True
In [28]:
arma_t.isstationary
Out[28]:
False
  • What does this mean?
In [29]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(arma_t.generate_sample(nsample=50));
In [30]:
arparams = np.array([1, .35, -.15, .55, .1])
maparams = np.array([1, .65])
arma_t = ArmaProcess(arparams, maparams)
arma_t.isstationary
Out[30]:
True
In [31]:
arma_rvs = arma_t.generate_sample(nsample=500, burnin=250, scale=2.5)
In [32]:
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(arma_rvs, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(arma_rvs, lags=40, ax=ax2)
  • For mixed ARMA processes the Autocorrelation function is a mixture of exponentials and damped sine waves after (q-p) lags.
  • The partial autocorrelation function is a mixture of exponentials and dampened sine waves after (p-q) lags.
In [33]:
arma11 = sm.tsa.ARMA(arma_rvs, (1,1)).fit()
resid = arma11.resid
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'))
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:650: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  elif issubdtype(paramsdtype, complex):
            AC           Q      Prob(>Q)
lag                                     
1.0   0.254921   32.687678  1.082211e-08
2.0  -0.172416   47.670747  4.450706e-11
3.0  -0.420945  137.159392  1.548466e-29
4.0  -0.046875  138.271301  6.617702e-29
5.0   0.103240  143.675908  2.958721e-29
6.0   0.214864  167.132999  1.823719e-33
7.0  -0.000889  167.133401  1.009206e-32
8.0  -0.045418  168.185752  3.094835e-32
9.0  -0.061445  170.115802  5.837215e-32
10.0  0.034623  170.729855  1.958737e-31
11.0  0.006351  170.750556  8.267054e-31
12.0 -0.012882  170.835909  3.220233e-30
13.0 -0.053959  172.336547  6.181196e-30
14.0 -0.016606  172.478965  2.160215e-29
15.0  0.051742  173.864487  4.089546e-29
16.0  0.078917  177.094281  3.217936e-29
17.0 -0.001834  177.096029  1.093168e-28
18.0 -0.101604  182.471937  3.103823e-29
19.0 -0.057342  184.187772  4.624066e-29
20.0  0.026975  184.568286  1.235670e-28
21.0  0.062359  186.605963  1.530259e-28
22.0 -0.009400  186.652365  4.548194e-28
23.0 -0.068037  189.088185  4.562011e-28
24.0 -0.035566  189.755202  9.901094e-28
25.0  0.095679  194.592622  3.354289e-28
26.0  0.065650  196.874877  3.487622e-28
27.0 -0.018404  197.054613  9.008747e-28
28.0 -0.079244  200.394009  5.773712e-28
29.0  0.008499  200.432501  1.541386e-27
30.0  0.053372  201.953775  2.133192e-27
31.0  0.074816  204.949395  1.550161e-27
32.0 -0.071187  207.667242  1.262288e-27
33.0 -0.088145  211.843156  5.480815e-28
34.0 -0.025283  212.187450  1.215228e-27
35.0  0.125690  220.714898  8.231608e-29
36.0  0.142724  231.734119  1.923081e-30
37.0  0.095768  236.706161  5.937779e-31
38.0 -0.084744  240.607804  2.890885e-31
39.0 -0.150126  252.878984  3.962998e-33
40.0 -0.083767  256.707741  1.996171e-33
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:577: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
In [34]:
arma41 = sm.tsa.ARMA(arma_rvs, (4,1)).fit()
resid = arma41.resid
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'))
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:646: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:650: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  elif issubdtype(paramsdtype, complex):
            AC          Q  Prob(>Q)
lag                                
1.0  -0.007888   0.031301  0.859570
2.0   0.004132   0.039906  0.980245
3.0   0.018103   0.205415  0.976710
4.0  -0.006760   0.228539  0.993948
5.0   0.018120   0.395025  0.995466
6.0   0.050688   1.700448  0.945087
7.0   0.010252   1.753955  0.972196
8.0  -0.011206   1.818018  0.986092
9.0   0.020292   2.028518  0.991009
10.0  0.001029   2.029060  0.996113
11.0 -0.014035   2.130168  0.997984
12.0 -0.023858   2.422925  0.998427
13.0 -0.002108   2.425216  0.999339
14.0 -0.018783   2.607429  0.999590
15.0  0.011316   2.673698  0.999805
16.0  0.042159   3.595419  0.999443
17.0  0.007943   3.628205  0.999734
18.0 -0.074311   6.503855  0.993686
19.0 -0.023379   6.789067  0.995256
20.0  0.002398   6.792073  0.997313
21.0  0.000487   6.792198  0.998516
22.0  0.017952   6.961435  0.999024
23.0 -0.038576   7.744466  0.998744
24.0 -0.029816   8.213249  0.998859
25.0  0.077850  11.415824  0.990675
26.0  0.040408  12.280448  0.989479
27.0 -0.018612  12.464276  0.992262
28.0 -0.014764  12.580187  0.994586
29.0  0.017650  12.746191  0.996111
30.0 -0.005486  12.762264  0.997504
31.0  0.058256  14.578545  0.994614
32.0 -0.040840  15.473084  0.993887
33.0 -0.019493  15.677309  0.995393
34.0  0.037269  16.425466  0.995214
35.0  0.086212  20.437448  0.976296
36.0  0.041271  21.358847  0.974774
37.0  0.078704  24.716878  0.938948
38.0 -0.029729  25.197056  0.944895
39.0 -0.078397  28.543388  0.891179
40.0 -0.014466  28.657579  0.909268
/ws/builds/jenkins/ws/du3/components/statsmodels/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/tsa/kalmanf/kalmanfilter.py:577: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if issubdtype(paramsdtype, float):

Exercise: How good of in-sample prediction can you do for another series, say, CPI

In [35]:
macrodta = sm.datasets.macrodata.load_pandas().data
macrodta.index = pd.Index(sm.tsa.datetools.dates_from_range('1959Q1', '2009Q3'))
cpi = macrodta["cpi"]

Hint:

In [36]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax = cpi.plot(ax=ax);
ax.legend();

P-value of the unit-root test, resoundly rejects the null of no unit-root.

In [37]:
print(sm.tsa.adfuller(cpi)[1])
0.990432818833742