0001-01-01 2045 words 5 minutes

Contents

时间序列

概述

分类
- 研究对象：
  - 一元时间序列
  - 多元时间序列
- 时间连续性
  - 离散时间序列
  - 连续时间序列
- 序列的特性
  - 平稳时间序列
  - 非平稳时间序列
分析方法
- 效应分解法
  - 将时间序列分解为趋势、周期、随机三个部分，并对前两个部分使用曲线进行拟合
- 平稳时间序列模型
  - 通过分析前后观测点之间的相关关系构建动态微分方程，用于预测

效应分解法

加法模型

$$ 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()