text: # coding: utf-8
  
  # DO NOT EDIT
  # Autogenerated from the notebook tsa_arma_0.ipynb.
  # Edit the notebook and then sync the output with this file.
  #
  # flake8: noqa
  # DO NOT EDIT
  
  # # 自回归移动平均线（ARMA）：太阳黑子数据
  
  import numpy as np
  from scipy import stats
  import pandas as pd
  import matplotlib.pyplot as plt
  
  import statsmodels.api as sm
  
  from statsmodels.graphics.api import qqplot
  
  # ## 太阳黑子数据
  
  print(sm.datasets.sunspots.NOTE)
  
  dta = sm.datasets.sunspots.load_pandas().data
  
  dta.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(disp=False)
  print(arma_mod20.params)
  
  arma_mod30 = sm.tsa.ARMA(dta, (3, 0)).fit(disp=False)
  
  print(arma_mod20.aic, arma_mod20.bic, arma_mod20.hqic)
  
  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.resid
  
  stats.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.loc['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)
  
  # ### 练习: 您能否更好的拟合 Sunspots 模型？ （提示：sm.tsa.AR 有一种方法 select_order）
  
  # ### 模拟ARMA（4,1）：模型识别困难
  
  from statsmodels.tsa.arima_process import ArmaProcess
  
  np.random.seed(1234)
  # 包括零滞后
  arparams = np.array([1, .75, -.65, -.55, .9])
  maparams = np.array([1, .65])
  
  # 让我们确保这个模型是可估计的。
  
  arma_t = ArmaProcess(arparams, maparams)
  
  arma_t.isinvertible
  
  arma_t.isstationary
  
  # * 这是什么意思？
  
  fig = plt.figure(figsize=(12, 8))
  ax = fig.add_subplot(111)
  ax.plot(arma_t.generate_sample(nsample=50))
  
  arparams = np.array([1, .35, -.15, .55, .1])
  maparams = np.array([1, .65])
  arma_t = ArmaProcess(arparams, maparams)
  arma_t.isstationary
  
  arma_rvs = arma_t.generate_sample(nsample=500, burnin=250, scale=2.5)
  
  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)
  
  # * 对于ARMA混合过程，自相关函数是（q-p）滞后后的指数和阻尼正弦波的混合。
  # * 部分自相关函数是（p-q）滞后后的指数和阻尼正弦波的混合。
  
  arma11 = sm.tsa.ARMA(arma_rvs, (1, 1)).fit(disp=False)
  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'))
  
  arma41 = sm.tsa.ARMA(arma_rvs, (4, 1)).fit(disp=False)
  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'))
  
  # ### 练习：您对另一个系列（例如CPI）的样本内预测有多好
  
  macrodta = sm.datasets.macrodata.load_pandas().data
  macrodta.index = pd.Index(
      sm.tsa.datetools.dates_from_range('1959Q1', '2009Q3'))
  cpi = macrodta["cpi"]
  
  # #### 提示:
  
  fig = plt.figure(figsize=(12, 8))
  ax = fig.add_subplot(111)
  ax = cpi.plot(ax=ax)
  ax.legend()
  
  # 单位根检验的 P-value，强烈地拒绝了一个单位根的 null值。
  
  print(sm.tsa.adfuller(cpi)[1])
  