ARMA模型的Python求解¶
试利用太阳黑子个数文件 sunspots.csv, 建立适 当的 ARMA 模型, 并预测 1989 年太阳黑子个数。 解 可以使用 statsmodels.api.tsa.ARMA()函数来拟合 ARMA 模型, 下面先初步地使用 ARMA $(9,1)$ 模型来拟合数据。得到 1989 年太阳黑子预测 值为 141 个。
代码
import pandas as pd, numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.rc('font',family='SimHei'); plt.rc('font',size=16)
d=pd.read_csv('sunspots.csv',usecols=['counts'])
md=sm.tsa.ARMA(d,(9,1)).fit()
years=np.arange(1700,1989) #已知观测值的年代
dhat=md.predict()
plt.plot(years[-20:],d.values[-20:],'o-k')
plt.plot(years[-20:],dhat.values[-20:],'P--')
plt.legend(('原始观测值','预测值')); plt.show()
dnext=md.predict(d.shape[0],d.shape[0])
print(dnext) #显示下一期的预测值
下面给出一个完整的建模步骤。 第一步：画出原始数据的折线图如图 $18.2$ 所示, 初步确定观测数据是 平稳的。画出序列的自相关图和偏相关图如图 。
第二步: 利用 AIC 和 BIC 准则, 确定选择 ARMA(4,2), 利用 Python 软 件，求得模型的计算结果，残差取值及分布 。 第三步：利用得到的模型, 得到 1989 年太阳黑子预测值为 139 个。原始数据及其预测值对比。
代码
import pandas as pd, numpy as np
import statsmodels.api as sm
import matplotlib.pylab as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.rc('axes',unicode_minus=False)
plt.rc('font',family='SimHei'); plt.rc('font',size=16)
d=pd.read_csv('sunspots.csv'); dd=d['counts']
years=d['year'].values.astype(int)
plt.plot(years,dd.values,'-*'); plt.figure()
ax1=plt.subplot(121); plot_acf(dd,ax=ax1,title='自相关')
ax2=plt.subplot(122); plot_pacf(dd,ax=ax2,title='偏自相关')
for i in range(1,6):
for j in range(1,6):
md=sm.tsa.ARMA(dd,(i,j)).fit()
print([i,j,md.aic,md.bic])
zmd=sm.tsa.ARMA(dd,(4,2)).fit()
print(zmd.summary()) #显示模型的所有信息
residuals = pd.DataFrame(zmd.resid)
fig, ax = plt.subplots(1,2)
residuals.plot(title="残差", ax=ax[0])
residuals.plot(kind='kde', title='密度', ax=ax[1])
plt.legend(''); plt.ylabel('')
dhat=zmd.predict(); plt.figure()
plt.plot(years[-20:],dd.values[-20:],'o-k')
plt.plot(years[-20:],dhat.values[-20:],'P--')
plt.legend(('原始观测值','预测值'))
dnext=zmd.predict(d.shape[0],d.shape[0])
print(dnext) #显示下一期的预测值
plt.show()
