关于GDP的ARIMA模型

这是一个时间序列课程论文项目
做了个arima
数据来源是国家统计局的gdp数据
可自行下载

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacf
import time

# 字体设置
font2 = {'family' : 'Times New Roman',
'weight' : 'normal',
'size'   : 20,
}

print(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())) 

首先，导入数据：

GDP_data=pd.read_excel('GDP年度数据.xls') 
# 切片出gdp数据
gdp_np=np.array(GDP_data.iloc[3,1:])[::-1]
# 创建时间序列
index=pd.date_range('1968','2021',freq='y')
gdp_serise=pd.Series(gdp_np,index=index)

做出GDP图像：

fig=plt.figure(figsize=(12,6))

plt.plot(gdp_serise)
plt.ylabel('GDP/(100 million yuan)',font2)
plt.xlabel('Year',font2)
plt.legend()
plt.locator_params(axis='x', nbins=12)
plt.tick_params(labelsize=16)
#plt.savefig('gdp.png', dpi=800, bbox_inches='tight')
plt.show()

对数据取ln，并创建时间序列：

gdp_np_log = []
for gdp in gdp_np:
    gdp_np_log.append(np.log(gdp))

gdp_series_log = pd.Series(data=np.array(gdp_np_log).ravel(),index=index)

做出ln(GDP)图像：

fig=plt.figure(figsize=(12,6))
plt.plot(gdp_series_log)
plt.ylabel('ln(GDP)',font2)
plt.xlabel('Year',font2)
plt.legend()
plt.locator_params(axis='x', nbins=12)
plt.tick_params(labelsize=16)
#plt.savefig('ln(gdp).png', dpi=800, bbox_inches='tight')
plt.show()

肯定不平稳，可以做adf检验，可以用eviews做，这里因为python做出来的和eviews有区别，所以未展示python做法，差分：

d_gdp_series_log=gdp_series_log.diff()

做出差分后图像：

fig=plt.figure(figsize=(12,6))
plt.plot(d_gdp_series_log['1969':'2019'])
plt.ylabel('diff_ln(GDP)',font2)
plt.xlabel('Year',font2)
plt.legend()
plt.locator_params(axis='x', nbins=12)
plt.tick_params(labelsize=16)
#plt.savefig('diff_ln(gdp).png', dpi=800, bbox_inches='tight')
plt.show()

画ACF、PACF：

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(d_gdp_series_log[1:],lags=24,ax=ax1)

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_pacf(d_gdp_series_log[1:],lags=24,ax=ax1)

ln(gdp)使用ARMA(1,0)，即相当于对GDP使用ARIMA(1,1,0):

model = ARIMA(gdp_series_log,order=(1,1,0))
# 拟合模型
result=model.fit()
print(result.summary()) 
# 预测结果
pred = result.predict()

由于pred为差分后的预测结果，故需要进行反向差分，再反向取ln：

# 反向差分
pred[0]=pred[0]+gdp_series_log[0]
for i in range(len(pred)-1):
    pred[i+1]=pred[i]+pred[i+1]
   
 # 反向ln
for i in range(len(pred)):
    pred[i]=np.exp(pred[i])

作图：

# 作图
fig=plt.figure(figsize=(12,6))
plt.plot(pred[1:],'r--',label='y_hat')
plt.plot(gdp_serise[1:],label='y')
plt.ylabel('GDP/(100 million yuan)',font2)
plt.xlabel('Year',font2)
plt.legend(fontsize='x-large')
plt.locator_params(axis='x', nbins=12)
plt.tick_params(labelsize=16)
#plt.savefig('arima预测.png', dpi=800, bbox_inches='tight')
plt.show()

定义预测精度评价指标，RMSE、MAE、MAPE：

from  sklearn import  metrics
def GetRMSE(y_hat,y_test):
    sum = np.sqrt(metrics.mean_squared_error(y_test, y_hat))
    return  sum

def GetMAE(y_hat,y_test):
    sum = metrics.mean_absolute_error(y_test, y_hat)
    return  sum

def GetMAPE(y_hat,y_test):
    sum = np.mean(np.abs((y_hat - y_test) / y_test)) * 100
    return sum

输出预测精度（样本内预测）：

print(GetRMSE(pred[1:-2],gdp_serise[2:]))
print(GetMAE(pred[1:-2],gdp_serise[2:]))
print(GetMAPE(pred[1:-2],gdp_serise[2:]))

还可以用eviews做一个残差分析图。

over！