Autoregressive Moving Average (ARMA): Sunspots data

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

import statsmodels.api as sm
[3]:
from statsmodels.graphics.api import qqplot

Sunspots Data

[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.

[5]:
dta = sm.datasets.sunspots.load_pandas().data
[6]:
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))
del dta["YEAR"]
[7]:
dta.plot(figsize=(12,8));
../../../_images/examples_notebooks_generated_tsa_arma_0_8_0.png
[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)
../../../_images/examples_notebooks_generated_tsa_arma_0_9_0.png
[9]:
arma_mod20 = sm.tsa.ARMA(dta, (2,0)).fit(disp=False)
print(arma_mod20.params)
const                49.659426
ar.L1.SUNACTIVITY     1.390656
ar.L2.SUNACTIVITY    -0.688571
dtype: float64
/usr/lib/python3/dist-packages/statsmodels/tsa/base/tsa_model.py:159: ValueWarning: No frequency information was provided, so inferred frequency A-DEC will be used.
  warnings.warn('No frequency information was'
[10]:
arma_mod30 = sm.tsa.ARMA(dta, (3,0)).fit(disp=False)
/usr/lib/python3/dist-packages/statsmodels/tsa/base/tsa_model.py:159: ValueWarning: No frequency information was provided, so inferred frequency A-DEC will be used.
  warnings.warn('No frequency information was'
[11]:
print(arma_mod20.aic, arma_mod20.bic, arma_mod20.hqic)
2622.6363380637877 2637.5697031713785 2628.6067259090337
[12]:
print(arma_mod30.params)
const                49.749962
ar.L1.SUNACTIVITY     1.300810
ar.L2.SUNACTIVITY    -0.508093
ar.L3.SUNACTIVITY    -0.129650
dtype: float64
[13]:
print(arma_mod30.aic, arma_mod30.bic, arma_mod30.hqic)
2619.403628696482 2638.070335080971 2626.8666135030394
  • Does our model obey the theory?

[14]:
sm.stats.durbin_watson(arma_mod30.resid.values)
[14]:
1.9564808699288254
[15]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax = arma_mod30.resid.plot(ax=ax);
../../../_images/examples_notebooks_generated_tsa_arma_0_17_0.png
[16]:
resid = arma_mod30.resid
[17]:
stats.normaltest(resid)
[17]:
NormaltestResult(statistic=49.84501385314911, pvalue=1.5006961438664718e-11)
[18]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
../../../_images/examples_notebooks_generated_tsa_arma_0_20_0.png
[19]:
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)
../../../_images/examples_notebooks_generated_tsa_arma_0_21_0.png
[20]:
r,q,p = sm.tsa.acf(resid.values.squeeze(), fft=True, 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.009179   0.026286  8.712032e-01
2.0   0.041793   0.573039  7.508724e-01
3.0  -0.001335   0.573599  9.024488e-01
4.0   0.136089   6.408919  1.706205e-01
5.0   0.092468   9.111828  1.046860e-01
6.0   0.091948  11.793245  6.674343e-02
7.0   0.068748  13.297202  6.518981e-02
8.0  -0.015020  13.369230  9.976130e-02
9.0   0.187592  24.641907  3.393913e-03
10.0  0.213718  39.321990  2.229478e-05
11.0  0.201082  52.361130  2.344957e-07
12.0  0.117182  56.804179  8.574293e-08
13.0 -0.014055  56.868316  1.893910e-07
14.0  0.015398  56.945555  3.997674e-07
15.0 -0.024967  57.149309  7.741499e-07
16.0  0.080916  59.296762  6.872184e-07
17.0  0.041138  59.853731  1.110947e-06
18.0 -0.052021  60.747421  1.548436e-06
19.0  0.062496  62.041684  1.831648e-06
20.0 -0.010302  62.076971  3.381251e-06
21.0  0.074453  63.926644  3.193596e-06
22.0  0.124955  69.154761  8.978387e-07
23.0  0.093162  72.071023  5.799805e-07
24.0 -0.082152  74.346677  4.713033e-07
25.0  0.015695  74.430032  8.289070e-07
26.0 -0.025037  74.642891  1.367289e-06
27.0 -0.125861  80.041135  3.722581e-07
28.0  0.053225  81.009970  4.716296e-07
29.0 -0.038693  81.523795  6.916659e-07
30.0 -0.016904  81.622214  1.151665e-06
31.0 -0.019296  81.750927  1.868772e-06
32.0  0.104990  85.575052  8.927992e-07
33.0  0.040086  86.134554  1.247514e-06
34.0  0.008829  86.161797  2.047833e-06
35.0  0.014588  86.236434  3.263819e-06
36.0 -0.119329  91.248884  1.084458e-06
37.0 -0.036665  91.723852  1.521929e-06
38.0 -0.046193  92.480501  1.938742e-06
39.0 -0.017768  92.592869  2.990691e-06
40.0 -0.006220  92.606692  4.697001e-06
  • This indicates a lack of fit.

  • In-sample dynamic prediction. How good does our model do?

[21]:
predict_sunspots = arma_mod30.predict('1990', '2012', dynamic=True)
print(predict_sunspots)
1990-12-31    167.047431
1991-12-31    140.993034
1992-12-31     94.859162
1993-12-31     46.860956
1994-12-31     11.242636
1995-12-31     -4.721255
1996-12-31     -1.166890
1997-12-31     16.185700
1998-12-31     39.021885
1999-12-31     59.449876
2000-12-31     72.170156
2001-12-31     75.376808
2002-12-31     70.436491
2003-12-31     60.731624
2004-12-31     50.201834
2005-12-31     42.076060
2006-12-31     38.114314
2007-12-31     38.454663
2008-12-31     41.963831
2009-12-31     46.869301
2010-12-31     51.423276
2011-12-31     54.399737
2012-12-31     55.321713
Freq: A-DEC, dtype: float64
[22]:
fig, ax = plt.subplots(figsize=(12, 8))
ax = dta.loc['1950':].plot(ax=ax)
fig = arma_mod30.plot_predict('1990', '2012', dynamic=True, ax=ax, plot_insample=False)
../../../_images/examples_notebooks_generated_tsa_arma_0_26_0.png
[23]:
def mean_forecast_err(y, yhat):
    return y.sub(yhat).mean()
[24]:
mean_forecast_err(dta.SUNACTIVITY, predict_sunspots)
[24]:
5.636930696153281

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

[25]:
from statsmodels.tsa.arima_process import ArmaProcess
[26]:
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.

[27]:
arma_t = ArmaProcess(arparams, maparams)
[28]:
arma_t.isinvertible
[28]:
True
[29]:
arma_t.isstationary
[29]:
False
  • What does this mean?

[30]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(arma_t.generate_sample(nsample=50));
../../../_images/examples_notebooks_generated_tsa_arma_0_38_0.png
[31]:
arparams = np.array([1, .35, -.15, .55, .1])
maparams = np.array([1, .65])
arma_t = ArmaProcess(arparams, maparams)
arma_t.isstationary
[31]:
True
[32]:
arma_rvs = arma_t.generate_sample(nsample=500, burnin=250, scale=2.5)
[33]:
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)
../../../_images/examples_notebooks_generated_tsa_arma_0_41_0.png
  • 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.

[34]:
arma11 = sm.tsa.ARMA(arma_rvs, (1,1)).fit(disp=False)
resid = arma11.resid
r,q,p = sm.tsa.acf(resid, fft=True, 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.254921   32.687694  1.082202e-08
2.0  -0.172416   47.670772  4.450649e-11
3.0  -0.420945  137.159409  1.548453e-29
4.0  -0.046875  138.271320  6.617642e-29
5.0   0.103240  143.675931  2.958688e-29
6.0   0.214864  167.133017  1.823703e-33
7.0  -0.000889  167.133419  1.009197e-32
8.0  -0.045418  168.185772  3.094806e-32
9.0  -0.061445  170.115821  5.837164e-32
10.0  0.034623  170.729873  1.958720e-31
11.0  0.006351  170.750574  8.266983e-31
12.0 -0.012882  170.835927  3.220205e-30
13.0 -0.053959  172.336565  6.181144e-30
14.0 -0.016606  172.478983  2.160197e-29
15.0  0.051742  173.864506  4.089511e-29
16.0  0.078917  177.094299  3.217908e-29
17.0 -0.001834  177.096047  1.093158e-28
18.0 -0.101604  182.471956  3.103796e-29
19.0 -0.057342  184.187791  4.624025e-29
20.0  0.026975  184.568306  1.235659e-28
21.0  0.062359  186.605982  1.530245e-28
22.0 -0.009400  186.652384  4.548155e-28
23.0 -0.068037  189.088205  4.561969e-28
24.0 -0.035566  189.755221  9.901006e-28
25.0  0.095679  194.592642  3.354261e-28
26.0  0.065650  196.874897  3.487591e-28
27.0 -0.018404  197.054634  9.008666e-28
28.0 -0.079244  200.394029  5.773663e-28
29.0  0.008499  200.432521  1.541373e-27
30.0  0.053372  201.953796  2.133173e-27
31.0  0.074816  204.949414  1.550148e-27
32.0 -0.071187  207.667261  1.262278e-27
33.0 -0.088145  211.843176  5.480766e-28
34.0 -0.025283  212.187470  1.215217e-27
35.0  0.125690  220.714917  8.231541e-29
36.0  0.142724  231.734141  1.923063e-30
37.0  0.095768  236.706182  5.937724e-31
38.0 -0.084744  240.607825  2.890859e-31
39.0 -0.150126  252.879008  3.962957e-33
40.0 -0.083767  256.707765  1.996150e-33
[35]:
arma41 = sm.tsa.ARMA(arma_rvs, (4,1)).fit(disp=False)
resid = arma41.resid
r,q,p = sm.tsa.acf(resid, fft=True, 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.007889   0.031302  0.859569
2.0   0.004132   0.039907  0.980244
3.0   0.018103   0.205418  0.976710
4.0  -0.006760   0.228541  0.993948
5.0   0.018120   0.395025  0.995466
6.0   0.050688   1.700445  0.945087
7.0   0.010252   1.753952  0.972197
8.0  -0.011206   1.818014  0.986092
9.0   0.020292   2.028515  0.991009
10.0  0.001029   2.029058  0.996113
11.0 -0.014035   2.130166  0.997984
12.0 -0.023858   2.422923  0.998427
13.0 -0.002108   2.425214  0.999339
14.0 -0.018783   2.607427  0.999590
15.0  0.011316   2.673696  0.999805
16.0  0.042159   3.595416  0.999443
17.0  0.007943   3.628201  0.999734
18.0 -0.074311   6.503850  0.993686
19.0 -0.023379   6.789062  0.995256
20.0  0.002398   6.792069  0.997313
21.0  0.000487   6.792193  0.998516
22.0  0.017952   6.961430  0.999024
23.0 -0.038576   7.744460  0.998744
24.0 -0.029816   8.213243  0.998859
25.0  0.077850  11.415817  0.990675
26.0  0.040408  12.280441  0.989479
27.0 -0.018612  12.464268  0.992262
28.0 -0.014764  12.580179  0.994586
29.0  0.017650  12.746183  0.996111
30.0 -0.005486  12.762256  0.997504
31.0  0.058256  14.578537  0.994614
32.0 -0.040840  15.473076  0.993887
33.0 -0.019493  15.677299  0.995393
34.0  0.037269  16.425456  0.995214
35.0  0.086212  20.437440  0.976296
36.0  0.041271  21.358837  0.974774
37.0  0.078704  24.716868  0.938949
38.0 -0.029729  25.197044  0.944895
39.0 -0.078397  28.543372  0.891179
40.0 -0.014466  28.657562  0.909269

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

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

Hint:

[37]:
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax = cpi.plot(ax=ax);
ax.legend();
../../../_images/examples_notebooks_generated_tsa_arma_0_48_0.png

P-value of the unit-root test, resoundingly rejects the null of a unit-root.

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