#!/usr/bin/env pythonfrom __future__ import print_functionimport numpy as npfrom scipy import statsimport pandas as pdimport matplotlibmatplotlib.style.use('ggplot')import matplotlib.pyplot as pltimport statsmodels.api as smfrom statsmodels.graphics.api import qqplotprint(sm.datasets.sunspots.NOTE)dta = sm.datasets.sunspots.load_pandas().datadta.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))del dta["YEAR"]dta.plot(figsize=(12,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)arma_mod20 = sm.tsa.ARMA(dta, (2,0)).fit()print(arma_mod20.params)print(arma_mod20.aic, arma_mod20.bic, arma_mod20.hqic)arma_mod30 = sm.tsa.ARMA(dta, (3,0)).fit()print(arma_mod30.params)print(arma_mod30.aic, arma_mod30.bic, arma_mod30.hqic)sm.stats.durbin_watson(arma_mod30.resid.values)fig = plt.figure(figsize=(12,8))ax = fig.add_subplot(111)ax = arma_mod30.resid.plot(ax=ax);  ###一定记得加分号，要不然后面的图片显示不了resid = arma_mod30.residstats.normaltest(resid)fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(111)fig = qqplot(resid, line='q', ax=ax, fit=True)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)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'))predict_sunspots = arma_mod30.predict('1990', '2012', dynamic=True)print(predict_sunspots)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)def mean_forecast_err(y, yhat):    return y.sub(yhat).mean()mean_forecast_err(dta.SUNACTIVITY, predict_sunspots)##############""" ADF test """import numpy as npimport statsmodels.tsa.stattools as tsx = np.array([1,2,3,4,3,4,2,3])result = ts.adfuller(x, 1)print("Test statistic", result[0])print("p-value", result[1])print("Lags Used", result[2])print("Number of observations Used", result[3])print("Critical Value(1%)", result[4]['1%'])print("Critical Value(5%)", result[4]['5%'])print("Critical Value(10%)", result[4]['10%'])