时间序列


import pandas as pd, numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf
from matplotlib.ticker import MaxNLocator
def acf_pacf(series, lags=40, title=None, figsize=(10,6)):
    #ACF and PACF plots
    #
求自相关和偏自相关函数
    # series
输入的时间序列
    # lags
自相关和偏自相关函数的滞后取值范围
    fig = plt.figure(figsize=figsize)
    ax1 = fig.add_subplot(211)
    ax1.set_xlabel('lags')
    ax1.set_ylabel('Autocorrelation')
    # x坐标为整数
    ax1.xaxis.set_major_locator(MaxNLocator(integer=True))
    plot_acf(series.dropna(), lags=lags, ax=ax1, title='ACF of %s'%title)
    ax2 = fig.add_subplot(212)
    ax2.set_xlabel('lags')
    ax2.set_ylabel('Partial Autocorrelation')
    ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
    plot_pacf(series.dropna(), lags=lags, ax=ax2, title='PACF of %s'%title)
    fig.tight_layout()
    
def test_stationarity(timeseries, window=10):
    #
求移动平均和移动标准差
    # window为选取的时间窗口个数
    rolmean = timeseries.rolling(window=window, center=False).mean()
    rolstd = timeseries.rolling(window=window, center=False).std()
 
    #
将以周重取样后的时间序列、移动平均和移动标准差制成折线图
#    orig = plt.plot(timeseries, color='blue', label='Original')
#    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
#    std = plt.plot(rolstd, color='black', label='Rolling Std')
#    plt.legend(loc='best')
#    plt.title('Rolling Mean & Standard Deviation')
#    plt.show(block=False)
 
    #
用Augmented Dickey-Fuller检验测试时间序列稳定性:
    print('Results of Augmented Dickey-Fuller Test:')
    #
使用减小AIC的办法估算ADF测试所需的滞后数
    dftest = adfuller(timeseries, autolag='AIC')
    #
将ADF测试结果、显著性概率、所用的滞后数和所用的观测数打印出来
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value',

                                             'Num Lags Used',

                                             'Num Observations Used'])
    for key, value in dftest[4].items():
        dfoutput['Critical Value (%s)' % key] = value
#    print(dfoutput)    
    
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.stats.diagnostic import acorr_ljungbox
def decompose_plot(series, title=''):
    decomposition = seasonal_decompose(series)
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    fig = decomposition.plot()
    fig.set_size_inches(10,6)
    fig.suptitle(title)
    fig.tight_layout()
    fig2 = acf_pacf(residual, title='Residuals', figsize=(10,4))
 
#
读取CSV
df = pd.read_csv('D:/jingjihu/portland-oregon-average-monthly-.csv')
#print(df.head(), '\nindex type:\n', type(df.index))
 
#
按'%Y-%m'格式转换CSV的日期,
为了方便处理,
将时间段转为timestamp
df['Month'] = pd.to_datetime(df['Month'], format='%Y-%m')
 
#
索引并resample为月
indexed_df = df.set_index('Month')
ts = indexed_df['riders']
ts = ts.resample('M').sum()
print(ts.head(), '\nindex type:\n', type(ts.index))
 
ts.plot(title='Monthly Num. of Ridership')
test_stationarity(ts)
decomposition = seasonal_decompose(ts)  
fig = plt.figure()  
fig = decomposition.plot()  
fig.set_size_inches(10, 6)
acf_pacf(decomposition.resid, title='Residuals')
ts_sdiff = ts - ts.shift(12)
print(ts_sdiff.shape)
test_stationarity(ts_sdiff.dropna())
 
decomposition = seasonal_decompose(ts_sdiff.dropna(), freq=12)  
fig = plt.figure()  
fig = decomposition.plot()  
fig.set_size_inches(10, 6)
acf_pacf(decomposition.resid, title='Residuals')
ts_diff_of_sdiff = ts_sdiff - ts_sdiff.shift(1)
test_stationarity(ts_diff_of_sdiff.dropna())
decomposition = seasonal_decompose(ts_diff_of_sdiff.dropna(), freq=12)  
fig = plt.figure()  
fig = decomposition.plot()  
fig.set_size_inches(10, 6)
acf_pacf(decomposition.resid, title='Residuals')
ts_log = np.log(ts)
ts_log_sdiff = ts_log - ts_log.shift(12)
ts_diff_of_log_sdiff = ts_log_sdiff - ts_log_sdiff.shift(1)
test_stationarity(ts_diff_of_log_sdiff.dropna())
decomposition = seasonal_decompose(ts_diff_of_log_sdiff.dropna(), freq=12)  
fig = plt.figure()  
fig = decomposition.plot()  
fig.set_size_inches(10, 6)
acf_pacf(decomposition.resid, title='Residuals')
from statsmodels.tsa.statespace.sarimax import SARIMAX
# SARIMA的参数: order = p, d, q, seasonal_order=P, D, Q, season的周期=12
model = SARIMAX(ts, order=(0,1,0),seasonal_order=(0,1,1,12))
#
已拟合的模型
fitted = model.fit()  
# ARIMA的结果
print(fitted.summary())
 
 比较拟合的数据和原始数据,
#ts.plot()
fitted.fittedvalues.plot(color='red')
 
#
拟合的数据和原始数据的残差
residuals = pd.DataFrame(fitted.resid)a
 
acf_pacf(residuals, title='Residuals', figsize=(10,4))
 
#
计算拟残差平方和(RSS)
plt.title('RSS: %.4f'% sum((fitted.fittedvalues-ts)**2))
 
#
残差的核密度估计
residuals.plot(kind='kde')
print('\n\nResiduals Summary:\n', residuals.describe())
#
plt.show()