# -*- coding:utf-8 -*-
import numpy as np
import pandas as pdfrom datetime import datetimeimport matplotlib.pylab as plt

# 读取数据，pd.read_csv默认生成DataFrame对象，需将其转换成Series对象df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date')df.index = pd.to_datetime(df.index)  # 将字符串索引转换成时间索引ts = df['x']  # 生成pd.Series对象# 查看数据格式ts.head()ts.head().index 

ts['1949-01-01']
ts[datetime(1949,1,1)]

ts['1949']

ts['1949-1' : '1949-6']

# -*- coding:utf-8 -*-
from statsmodels.tsa.stattools import adfuller
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 移动平均图
def draw_trend(timeSeries, size):f = plt.figure(facecolor='white')# 对size个数据进行移动平均rol_mean = timeSeries.rolling(window=size).mean()# 对size个数据进行加权移动平均rol_weighted_mean = pd.ewma(timeSeries, span=size)timeSeries.plot(color='blue', label='Original')rolmean.plot(color='red', label='Rolling Mean')rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')plt.legend(loc='best')plt.title('Rolling Mean')plt.show()def draw_ts(timeSeries):    f = plt.figure(facecolor='white')timeSeries.plot(color='blue')plt.show()'''　　Unit Root TestThe null hypothesis of the Augmented Dickey-Fuller is that there is a unitroot, with the alternative that there is no unit root. That is to say thebigger the p-value the more reason we assert that there is a unit root
'''
def testStationarity(ts):dftest = adfuller(ts)# 对上述函数求得的值进行语义描述dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])for key,value in dftest[4].items():dfoutput['Critical Value (%s)'%key] = valuereturn dfoutput# 自相关和偏相关图，默认阶数为31阶
def draw_acf_pacf(ts, lags=31):f = plt.figure(facecolor='white')ax1 = f.add_subplot(211)plot_acf(ts, lags=31, ax=ax1)ax2 = f.add_subplot(212)plot_pacf(ts, lags=31, ax=ax2)plt.show()

ts_log = np.log(ts)
test_stationarity.draw_ts(ts_log)

test_stationarity.draw_trend(ts_log, 12)

diff_12 = ts_log.diff(12)
diff_12.dropna(inplace=True)
diff_12_1 = diff_12.diff(1)
diff_12_1.dropna(inplace=True)
test_stationarity.testStationarity(diff_12_1)

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log, model="additive")trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

rol_mean = ts_log.rolling(window=12).mean()
rol_mean.dropna(inplace=True)
ts_diff_1 = rol_mean.diff(1)
ts_diff_1.dropna(inplace=True)
test_stationarity.testStationarity(ts_diff_1)

ts_diff_2 = ts_diff_1.diff(1)
ts_diff_2.dropna(inplace=True)

from statsmodels.tsa.arima_model import ARMA
model = ARMA(ts_diff_2, order=(1, 1))
result_arma = model.fit( disp=-1, method='css')

predict_ts = result_arma.predict()
# 一阶差分还原diff_shift_ts = ts_diff_1.shift(1)diff_recover_1 = predict_ts.add(diff_shift_ts)# 再次一阶差分还原
rol_shift_ts = rol_mean.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts)
# 移动平均还原
rol_sum = ts_log.rolling(window=11).sum()
rol_recover = diff_recover*12 - rol_sum.shift(1)
# 对数还原
log_recover = np.exp(rol_recover)
log_recover.dropna(inplace=True)

ts = ts[log_recover.index]  # 过滤没有预测的记录plt.figure(facecolor='white')
log_recover.plot(color='blue', label='Predict')
ts.plot(color='red', label='Original')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))
plt.show()

# 差分操作
def diff_ts(ts, d):global shift_ts_list#  动态预测第二日的值时所需要的差分序列global last_data_shift_listshift_ts_list = []last_data_shift_list = []tmp_ts = tsfor i in d:last_data_shift_list.append(tmp_ts[-i])print last_data_shift_listshift_ts = tmp_ts.shift(i)shift_ts_list.append(shift_ts)tmp_ts = tmp_ts - shift_tstmp_ts.dropna(inplace=True)return tmp_ts# 还原操作
def predict_diff_recover(predict_value, d):if isinstance(predict_value, float):tmp_data = predict_valuefor i in range(len(d)):tmp_data = tmp_data + last_data_shift_list[-i-1]elif isinstance(predict_value, np.ndarray):tmp_data = predict_value[0]for i in range(len(d)):tmp_data = tmp_data + last_data_shift_list[-i-1]else:tmp_data = predict_valuefor i in range(len(d)):try:tmp_data = tmp_data.add(shift_ts_list[-i-1])except:raise ValueError('What you input is not pd.Series type!')tmp_data.dropna(inplace=True)return tmp_data

diffed_ts = diff_ts(ts_log, d=[12, 1])
model = arima_model(diffed_ts)
model.certain_model(1, 1)
predict_ts = model.properModel.predict()
diff_recover_ts = predict_diff_recover(predict_ts, d=[12, 1])
log_recover = np.exp(diff_recover_ts)

def proper_model(data_ts, maxLag):init_bic = sys.maxintinit_p = 0init_q = 0init_properModel = Nonefor p in np.arange(maxLag):for q in np.arange(maxLag):model = ARMA(data_ts, order=(p, q))try:results_ARMA = model.fit(disp=-1, method='css')except:continuebic = results_ARMA.bicif bic < init_bic:init_p = pinit_q = qinit_properModel = results_ARMAinit_bic = bicreturn init_bic, init_p, init_q, init_properModel

from dateutil.relativedelta import relativedeltadef _add_new_data(ts, dat, type='day'):

if type == 'day':new_index = ts.index[-1] + relativedelta(days=1)elif type == 'month':new_index = ts.index[-1] + relativedelta(months=1)ts[new_index] = datdef add_today_data(model, ts,  data, d, type='day'):_add_new_data(ts, data, type)  # 为原始序列添加数据# 为滞后序列添加新值d_ts = diff_ts(ts, d)model.add_today_data(d_ts[-1], type)def forecast_next_day_data(model, type='day'):if model == None:raise ValueError('No model fit before')fc = model.forecast_next_day_value(type)return predict_diff_recover(fc, [12, 1])

ts_train = ts_log[:'1956-12']
ts_test = ts_log['1957-1':]diffed_ts = diff_ts(ts_train, [12, 1])
forecast_list = []for i, dta in enumerate(ts_test):if i%7 == 0:model = arima_model(diffed_ts)model.certain_model(1, 1)forecast_data = forecast_next_day_data(model, type='month')forecast_list.append(forecast_data)add_today_data(model, ts_train, dta, [12, 1], type='month')predict_ts = pd.Series(data=forecast_list, index=ts['1957-1':].index)log_recover = np.exp(predict_ts)original_ts = ts['1957-1':]