ARIMA, ACF, PACF
文章目录
偏自相关系数(Partial Auto-Correlation Coefficient, PAC) 和 自相关系数(Auto-Correlation Coefficient, ACF) 是两种常用的时间序列分析工具,可用于观察数据中是否存在某种模式或结构。如果数据中有明显的周期性变化,那么通过观察PAC和ACF,可以从这些指标中推断出数据的周期性特点。以下是一些常用的步骤:
- 首先要对原始数据进行可视化,以便大致了解其趋势。这可以通过折线图或柱状图等方式实现。
- 然后使用上述提到的sm.tsa.stattools.acf 和 sm.tsa.stattools.pacf 方法计算ACF和PACF,并绘制出对应的图示。注意,PACF反映了滞后项对当前项的影响,而忽略了之前的滞后项;而ACF则反映了所有的滞后项对当前项的影响。
- 观察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
时间序列分析步骤
-
白噪声检验:如果是白噪声,就没有分析意义了
-
平稳非白噪声序列:假定均值和方差是常数:AR MA ARMA
-
非平稳序列 : 差分 ARIMA SARIMA
-
单变量:ARMA>>GARCH
-
多变量:VAR>>MGARCH
平稳性检验 >> 白噪声检验 >> ACF、PACF >> 模型选择
突变点分析
CUSUM算法
C+ 和 C- 分别标记为上累积和和下累积和。如果超过阈值H,则判定为断点或者突变点
C
i
+
=
m
a
x
(
0
,
X
i
−
T
a
r
g
e
t
σ
+
C
i
−
1
+
−
K
)
C_i^+ = max(0,\frac{X_{i}-Target}{\sigma}+C_{i-1}^+-K)
Ci+=max(0,σXi−Target+Ci−1+−K)
C
i
−
=
m
i
n
(
0
,
X
i
−
T
a
r
g
e
t
σ
−
C
i
−
1
−
+
K
)
C_i^- = min(0,\frac{X_i-Target}{\sigma} - C_{i-1}^-+K)
Ci−=min(0,σXi−Target−Ci−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