时间序列
概述
-
分类
- 研究对象:
- 一元时间序列
- 多元时间序列
- 时间连续性
- 离散时间序列
- 连续时间序列
- 序列的特性
- 平稳时间序列
- 非平稳时间序列
- 研究对象:
-
分析方法
-
效应分解法
-
将时间序列分解为趋势、周期、随机三个部分,并对前两个部分使用曲线进行拟合
-
-
平稳时间序列模型
- 通过分析前后观测点之间的相关关系构建动态微分方程,用于预测
-
效应分解法
加法模型
$$ x_t = T_t + S_t + I_t $$
其中,$T$代表趋势效应,$S$代表季节效应,$I$代表随机效应
乘法模型
$$ x_t = T_t \times S_t \times I_t $$
代码实现
加法模型
import pandas as pd
from fbprophet import Prophet
import matplotlib.pyplot as plt
# 数据读取
df = pd.read_csv('AirPassengers.csv')
df['DATE'] = pd.to_datetime(df['DATE'])
###################################################
# Prophet 只能处理单变量的时间序列
# 导入的数据第一列必须是日期类型的变量,变量名必须为'ds'
# 导入的数据第二列必须是数值类型的变量,变量名必须为'y'
####################################################
df = df.rename(columns={'DATE':'ds',
'AIR':'y'})
# 模型训练
# 长期趋势部分的拟合函数为linear,可选项为 linear 和 logistic
# 置信区间为0.95 默认为0.8
model = Prophet(growth='linear',interval_width=0.95)
model.fit(df)
# 模型预测
# 'MS' 表示月度数据 即36个月
future_dates = model.make_future_dataframe(period=36,freq='MS')
forecast = model.predict(future_dates)
# 数据展示
forecast[['ds','yhat','yhat_lower','y_hat_upper']]
# uncertainty 表示是否绘制置信区间
model.plot(forecast,uncertainty=True)
乘法模型
如果需要做乘法模型,只需要对时间序列取自然对数
import pandas as pd
from fbprophet import Prophet
import matplotlib.pyplot as plt
# 数据读取
df = pd.read_csv('AirPassengers.csv')
df['DATE'] = pd.to_datetime(df['DATE'])
df['AIR'] = np.log(df['AIR'])
###################################################
# Prophet 只能处理单变量的时间序列
# 导入的数据第一列必须是日期类型的变量,变量名必须为'ds'
# 导入的数据第二列必须是数值类型的变量,变量名必须为'y'
####################################################
df = df.rename(columns={'DATE':'ds',
'AIR':'y'})
# 模型训练
# 长期趋势部分的拟合函数为linear,可选项为 linear 和 logistic
# 置信区间为0.95 默认为0.8
model = Prophet(growth='linear',interval_width=0.95)
model.fit(df)
# 模型预测
# 'MS' 表示月度数据 即36个月
future_dates = model.make_future_dataframe(period=36,freq='MS')
forecast = model.predict(future_dates)
forecast[['yhat','yhat_lower','y_hat_upper']] = np.exp(forecast[['yhat','yhat_lower','y_hat_upper']])
# 数据展示
forecast[['yhat','yhat_lower','y_hat_upper']]
# uncertainty 表示是否绘制置信区间
model.plot(forecast,uncertainty=True)
平稳时间序列模型
- 序列分类
- 严平稳时间序列
- 一个时间序列中,各期数据的联合概率分布与时间t无关
- 宽平稳时间序列
- 对任意时间,序列的均值、方差存在并为常数,且自协方差函数与自相关系数只与时间间隔k有关
- 严平稳时间序列
==平稳时间序列分析在于挖掘时间序列之间的关系,当时间序列中的关系被提取出来后,剩下的序列就应该是一个白噪声序列==
-
模型分类
- 自回归模型(Auto Regression Model, AR)
- 移动平均模型(Moving Average Model, MA)
- 自回归移动平均模型(Auto Regression Moving Average Model, ARMA)
-
自相关函数(Autocorrelation Function, ACF)
- 描述时间序列任意两个时间间隔k的相关系数
$$ ACF(k)=\rho_k=\frac{Cov(y_t,y_{t-k})}{Var(y_t)} $$
- 偏自相关函数(Partial Autocorrelation Function, PACF)
- 描述时间序列中在任意两个时间间隔k的时刻,去除1至k-1这个时间段中的其他数据的相关系数
$$ \rho^_k = Corr[y_t -E^(y_t|y_{t-1},\cdots,y_{t-k+1}),y_{t-k}] $$
AR 模型
AR模型认为时间序列档期观测值与前n期有线性关系,而与前n+1期无线性关系 $$ X_t = \alpha_o+\alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_nX_{t-n}+\epsilon_t\ \epsilon_t \sim N(0,\sigma^2) $$
MA 模型
MA模型认为,如果一个系统$t$时刻的响应$X_t$,与其以前时刻$t-1,t-2,\cdots$的响应$X_{t-1},X_{t-2},\cdots$无关,而与其以前时刻$t-1,t-2,\cdots,t-m$进入系统的扰动项$\epsilon_{t-1},\epsilon_{t-2},\cdots,\epsilon_{t-m}$存在一定的相关关系 $$ X_t = \mu+\beta_1X_{\epsilon-1}+\beta_2X_{\epsilon-2}+\cdots+\beta_mX_{\epsilon-m}+\epsilon_t $$
ARMA 模型
ARMA模型结合了AR模型和MA模型的共同特点,其认为序列是受到前期观测数据与系统扰动的共同影响 $$ X_t = \alpha_o+\alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_nX_{t-n}+\epsilon_t -\beta_1X_{\epsilon-1}-\beta_2X_{\epsilon-2}-\cdots-\beta_mX_{\epsilon-m} $$
模型的定阶与识别
使用ACF与PACF)
自相关系数 (ACF) 与偏自相关系数 (PACF) 可以用于判断平稳时间序列数据适合哪一种模型和阶数
AR(p) 模型 | MA(q) 模型 | ARMA(p,q) 模型 | |
---|---|---|---|
ACF | 拖尾 | q阶截尾 | 拖尾 |
PACF | p阶截尾 | 拖尾 | 拖尾 |
使用AIC与BIC
使用ACF与PACF进行定阶时,只能精确AR模型与MA模型的阶数,而对于ARMA模型无法确定阶数
AIC或BIC准则特别适用于ARMA模型,也适用于AR模型与MA模型
代码实现
import statsmodels.api as sm
# 读取数据,探索数据的平稳性
ts_simu200 = pd.read_csv('ts_simu200.csv',index_col='t')
datas = pd.date_range(start='2017/01/01',periods=200)
ts_simu200.set_index(datas,inplace=True)
dta = ts_simu200['AR1_a']
dta.plot(figsize(12,8))
# 定阶
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(dta,lags=20) #lags 表示滞后的阶数
fig = sm.graphics.tsa.plot_pacf(dta,lags=20)
plt.show()
# AR 建模
ar10 = sm.tsa.ARMA(dta,(1,0)).fit()
# 残差白噪声检验
resid = ar10.resid
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(resid,values.squeeze(),lags=20)
fig = sm.graphics.tsa.plot_pacf(resid,legs=20)
# AR模型预测
predict_dta = ar10.forecast(step=5)
import datatime
fig.ar10.plot_predict(pd.to_datetime('2017-01-01')+datetime.timedelta(days=190),
pd.to_datetime('2017-01-01')+datetime.timedelta(days=220),
dynamic=False,plot_insasmple=True)
plt.show()
非平稳时间序列模型
ARIMA 模型
很多时间序列模型数据都是非平稳的,直接以平稳时间序列分析方法进行分析是不合适的。可以通过差分等手段,将非平稳时间序列数据变成平稳时间序列,再采用ARIMA模型建模。 $$ \grad x_t^{(1)} = x_t-x_{t-1}\ \grad x_t^{(2)} = \grad x_t^{(1)} - \grad x_{t-1}^{(1)}\ \grad x_t^{(d)} = \grad x_t^{(d-1)} - \grad x_{t-1}^{(d-1)} $$
对于季节性的数据,可以采用一定周期s的差分运算(季节差分)提取季节信息 $$ \grad_sx_t = x_t-x_{t-s} $$
代码实现
一阶差分模型
from statsmodels.tsa.stattools import adfuller
import statsmodels.api as sm
# 读取数据,探索数据的平稳性
ts_simu200 = pd.read_csv('ts_simu200.csv',index_col='t')
datas = pd.date_range(start='2017/01/01',periods=200)
ts_simu200.set_index(datas,inplace=True)
dta = ts_simu200['ARIMA_110']
dta.plot(figsize(12,8))
result = adfuller(dta)
print(f'ADF Statistic:{result[0]}')
print(f'p-value:{result[1]}')
# p值为0.333,无法拒绝平稳假设,认为数据非平稳
# 对原数据进行差分
diff1 = dta.diff(1)
diff(1).plot(figsize(12,8))
# 定阶
ddta = diff1
ddta.dropna(inplace=True)
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(ddta,lags=20) #lags 表示滞后的阶数
fig = sm.graphics.tsa.plot_pacf(ddta,lags=20)
plt.show()
# ARIMA 建模
arima10 = sm.tsa.ARIMA(dta,(1,1,0)).fit()
# 残差白噪声检验
resid = arima10.resid
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(resid,values.squeeze(),lags=20)
fig = sm.graphics.tsa.plot_pacf(resid,legs=20)
# ARIMA模型预测
predict_dta = arima10.forecast(step=5)
import datatime
fig.arima10.plot_predict(pd.to_datetime('2017-01-01'),
pd.to_datetime('2017-01-01')+datetime.timedelta(days=220),
dynamic=False,plot_insasmple=True)
plt.show()
季节性差分模型
# 读取数据,探索数据的平稳性
sales_ts = pd.read_csv('tractor_sales.csv')
plt.figure(figsize=(10,5))
plt.plot(np.log(sales_ts))
sales_ts_log = np.log(sales_ts)
sales_ts_log.dropna(inplace=True)
sales_ts_log_diff = sales_ts_log.diff(1)
sales_ts_log_diff.dropna(inplace=True)
fig,axes = plt.subplots(1,2,sharey=False,sharex = False)
fig.set_figwidth(12)
fig.set_figheight(4)
smt.graphics.plot_acf(sales_ts_log_diff,lags=30,ax=axes[0],alpha=0.5)
smt.graphics.plot_pacf(sales_ts_log_diff,lags=30,ax=axes[0],alpha=0.5)
plt.tight_layout()
p = d = q = range(0,2)
pdq = list(itertools.product(p,d,q))
seasonal_pdq = [(x[0],x[1],x[2]) for x in pdq]
best_aic = np.inf
best_pdq = None
best_seasonal_pdq = None
temp_model = None
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
temp_model = sm.tsa.statespace.SARIMAX(sales_ts_log,
order = param,
seasonal_order = param_seasonal,
enforce_statopnarity = True,
enforce_invertibility = True)
results = temp_model.fit()
if results.aic < best_aic:
best_aic = results.aic
best_pdq = param
best_seasonal_pdq = param_seasonal
except:
continue
best_model = sm.tsa.statespace.SARIMAX(sales_ts_log,
order = (0,1,1),
seasonal_order = (1,0,1,12),
enforce_stationarity = True,
enforce_invertibility = True)
best_results = best_model.fit()
best_results.plot_diagnostics(lags=30,figsize=(16,32))
plt.show()
n_step = 36
pred_uc_95 = best_results.get_forecast(steps=nsteps,alpha=0.05)
pred_pr_95 = pred_uc_95.predicted_mean
pred_ci_95 = pred_uc_95.conf_int()
idx = pd.date_range(sales_ts.index[-1],periods=n_steps,freq='MS')
fc_95 = pd.DataFrame(np.column_stack([np.exp(pred_pr_95),np.exp(pred_ci_95)]),index=idx,
columns=['forecast','lower_ci_95','upper_ci_95'])
axis = sales_ts.plot(label='Observed',figsize=(15,6))
fc_95['forecast'].plot(ax=axis,label='Forecast',alpha=0.7)
axis.fill_between(fc_all.index,fc_95['lower_ci_95'],fc_95['upper_ci_95'],
color='k',alpha=0.25)
axis.set_xlabel('Years')
axis.set_ylabel('Tractor Sales')
plt.legend(loc='best')
plt.show()