时序数据分析：ARIMA, ACF, PACF

ARIMA, ACF, PACF

偏自相关系数(Partial Auto-Correlation Coefficient, PAC) 和 自相关系数(Auto-Correlation Coefficient, ACF) 是两种常用的时间序列分析工具，可用于观察数据中是否存在某种模式或结构。如果数据中有明显的周期性变化，那么通过观察PAC和ACF，可以从这些指标中推断出数据的周期性特点。以下是一些常用的步骤：

1. 首先要对原始数据进行可视化，以便大致了解其趋势。这可以通过折线图或柱状图等方式实现。
2. 然后使用上述提到的sm.tsa.stattools.acf 和 sm.tsa.stattools.pacf 方法计算ACF和PACF，并绘制出对应的图示。注意，PACF反映了滞后项对当前项的影响，而忽略了之前的滞后项；而ACF则反映了所有的滞后项对当前项的影响。
3. 观察PACF和ACF的图形形状和值，如果发现有较高的峰值，那么可能存在一定的周期性规律。一般来说，如果存在周期性，则需要关注以下几种现象： * 如果在ACF中发现有多个峰点（或至少有一个明显峰点），并且这些峰点呈周期性排列，那么说明数据可能存在一定的周期性； * 在PACF中也可能出现类似的峰值，但需要注意的是，如果PACF中的峰值在后面迅速衰减，则可能是随机过程的结果，而不是真正的周期性表现； * 同时也可以结合谱密度图等其他方法来确认周期性的存在。

ACF

从字面上的意思理解，自相关函数（Autocorrelation Function）就是用来计算时间序列自身的相关性的函数。
c o r r ( x t , x t − 1 ) = c o r r ( { x 2 , x 3 , . . . x n } , { x 1 , x 2 , . . . x n − 1 } ) corr(x_t,x_{t-1}) = corr(\{x_2,x_3,...x_n\},\{x_1,x_2,...x_{n-1}\}) corr(xt​,xt−1​)=corr({x2​,x3​,...xn​},{x1​,x2​,...xn−1​})

PACF

PACF，指偏自相关函数（Partial Autocorrelation Function）实际上计算的就是两个回归方程的残差序列的相关系数。

from sklearn.linear_model import LinearRegression

df['price-shift-1'] = df['price'].shift(1)
df['price-shift-2'] = df['price'].shift(2)

lr = LinearRegression()
lr.fit(df[['price-shift-1']].iloc[1:, :], df[['price']].iloc[1:, :])
res1 = df[['price']].iloc[1:, :] - lr.predict(df[['price-shift-1']].iloc[1:, :])  # 残差序列1

lr.fit(df[['price-shift-1']].iloc[2:, :], df[['price-shift-2']].iloc[2:, :])
res2 = df[['price-shift-2']].iloc[2:, :] - lr.predict(df[['price-shift-1']].iloc[2:, :])  # 残差序列2

print(np.corrcoef(res1['price'].values[1:], res2['price-shift-2'].values)[1, 0])

通过看图确定时间序列周期

import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf

# 计算ACF和PACF
acf = sm.tsa.stattools.acf(data_demo, nlags=30)
pacf = sm.tsa.stattools.pacf(data_demo, nlags=30)

#绘制ACF图和PACF图
rand=np.array([random.randrange(1,10)/10 for i in range(200)])
y=np.array([i for i in range(10)]*20)+rand

plot_acf(y,lags=100)
plot_pacf(y,lags=100)

傅里叶变换寻找周期

对于一条具备周期性的时间序列，它本身就很接近正弦波，所以它包含一个显著的正弦波，周期就是该正弦波的周期，而这个正弦波可以通过傅里叶变换找到，它将时序数据展开成三角函数的线性组合，得到每个展开项的系数，就是傅里叶系数。傅里叶系数越大，表明它所对应的正弦波的周期就越有可能是这份数据的周期。

先通过傅里叶变换找到可能的周期，再用自相关系数做排除，从而得到最可能的周期。

from scipy.fftpack import fft, fftfreq

fft_series = fft(data["value"].values)
power = np.abs(fft_series)
sample_freq = fftfreq(fft_series.size)

pos_mask = np.where(sample_freq > 0)
freqs = sample_freq[pos_mask]
powers = power[pos_mask]

top_k_seasons = 3
# top K=3 index
top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:]
top_k_power = powers[top_k_idxs]
fft_periods = (1 / freqs[top_k_idxs]).astype(int)

print(f"top_k_power: {top_k_power}")
print(f"fft_periods: {fft_periods}")

现在计算各个可能周期的自相关系数

from statsmodels.tsa.stattools import acf

# Expected time period
for lag in fft_periods:
    # lag = fft_periods[np.abs(fft_periods - time_lag).argmin()]
    acf_score = acf(data["value"].values, nlags=lag)[-1]
    print(f"lag: {lag} fft acf: {acf_score}")


expected_lags = np.array([timedelta(hours=12)/timedelta(minutes=5), timedelta(days=1)/timedelta(minutes=5), timedelta(days=7)/timedelta(minutes=5)]).astype(int)
for lag in expected_lags:
    acf_score = acf(data["value"].values, nlags=lag, fft=False)[-1]
    print(f"lag: {lag} expected acf: {acf_score}")

通过auto_arima确定合适的参数

from statsmodels.tsa.arima.model import ARIMA
from pmdarima import auto_arima

# 使用auto_arima自动寻找最佳ARIMA模型参数
model_auto = auto_arima(time_series, seasonal=False, m=1, trace=True)

# 拟合ARIMA模型
model = ARIMA(time_series, order=model_auto.order)
model_fit = model.fit()

# 进行预测
n_periods = 1  
forecast = model_fit.get_forecast(steps=n_periods).predicted_mean

对任意时间序列，当ACF图像呈现拖尾、且PACF图像呈现截尾状态时，当前时间序列适用AR模型，且PACF截尾的滞后阶数就是超参数p的理想值，如图：

对任意时间序列，当PACF图像呈现拖尾、且ACF图像呈现截尾状态时，当前时间序列适用MA模型，且ACF截尾的滞后阶数就是超参数q的理想值

https://zhuanlan.zhihu.com/p/634120397

时间序列分析步骤

1. 白噪声检验：如果是白噪声，就没有分析意义了

2. 平稳非白噪声序列：假定均值和方差是常数：AR MA ARMA

3. 非平稳序列 ： 差分 ARIMA SARIMA

4. 单变量：ARMA>>GARCH

5. 多变量：VAR>>MGARCH

平稳性检验 >> 白噪声检验 >> ACF、PACF >> 模型选择

突变点分析

CUSUM算法

C+ 和 C- 分别标记为上累积和和下累积和。如果超过阈值H，则判定为断点或者突变点
$$C_i^{+}=max(0,\frac{X_i-Target}{\sigma }+C_{i-1}^{+}-K)$$
$$C_i^{-}=min(0,\frac{X_i-Target}{\sigma }-C_{i-1}^{-}+K)$$

import pandas as pd  
import numpy as np  
  
  
def solve(k,h): 
	'''
	:param k: 参考值/允许偏差，防止对小波动过度敏感 
	:param h: 决策间隔/阈值，确定变点检测的灵敏度 :return: 变点的索引列表'''
    #读取数据  
    arr = pd.read_csv('finalDatasetWithRain_100.csv')['tmp']  
    #zscore标准化  
    arr = (arr - np.mean(arr)) / np.std(arr)  
    #初始化上累积和和下累积和  
    c1, c2 = [arr[0]], [arr[0]]  
    for i in arr[1:]:  
        c1.append(max(0, c1[-1] + i - k))  
        c2.append(min(0, c2[-1] - i + k))  
    c1, c2 = np.array(c1), np.array(c2)  
    #计算突变点的索引  
    bkp = np.where(np.where(c1 > h, True, False) + np.where(c2 < -h, True, False))  
    return bkp

EDM算法

import pandas as pd  
import numpy as np  
from scipy.spatial.distance import pdist, squareform  
  
def solve():  
    # 读取数据集  
    df = pd.read_csv('finalDatasetWithRain_100.csv')  
  
    # 提取海平面气压slp数据  
    slp = df['slp'].values.reshape(-1, 1)  
      
    # 构造距离矩阵  
    distance_matrix = squareform(pdist(slp, 'euclidean'))  
      
    # 使用 EDM 算法进行变点检测  
    max_diff = -np.inf  
      
    #初始化最大变点的位置列表  
    bkps =[0]  
    for i in range(1,tmp.shape[0]):  
        #变点一侧的距离中位数  
        median1 = np.median(distance_matrix[:i,:i])  
        #另个一变点的距离中位数  
        median2 = np.median(distance_matrix[i:,i:])  
        #变点两侧点之间的距离中位数  
        median_i = np.median(distance_matrix[:i,i:])  
          
        #两段之间的距离中位数和两段内部的距离中位数的差异  
        diff = median_i-(median1+median2)/2   
if  diff > max_diff:  
            max_diff = diff   
bkps = [i]  
        elif diff == max_diff:  
            bkps.append(i)  
              
    return bkps  
solve()

Pelt算法

import numpy as np
import ruptures as rpt
  
def solve(pen):
	#读取数据
	
	data=pd.read_csv('SBUX.US_M10_120.csv')[['close']]
	
	algo = rpt.Pelt().fit(data)
	
	bpks = algo.predict(pen=pen)
	
	return bpks