【时间序列模型】AR模型（原理剖析+MATLAB代码）

前言

时间序列分析方法包括频域分析方法和时域分析方法。时域分析方法是从序列自相关的角度揭示时间序列的发展规律，常用的模型如下：
自回归（AutoRegressive，AR）模型
移动平均（Moving Average，MA）模型
自回归移动平均（AutoRegressive Moving Average，ARMA）模型
差分自回归移动平均（AutoRegressive Integrated Moving Average，ARIMA）模型

一、AR模型的原理剖析

1.1 AR模型原理

对于一组观测数据的离散时间序列$X$={$x_1,x_2,...x_{N-1},x_N$}，每个时间序列的取值不仅与前面的各值相关，也与产生的干扰有关，因而该时间序列可用线性差分方程表示：
$$x_i={\phi }_1{x}_{i-1}+{\phi }_2{x}_{i-2}+...+{\phi }_n{x}_{i-n}+a_i-{\theta }_1{a}_{i-1}-{\theta }_2{a}_{i-2}-...-{\theta }_m{a}_{i-m}$$
式中，$i=1,2,...,N.$，$a_i为模型残差$。
当滑动平均系数${\theta }_j=0(j=1,2,...,m.)$时，即没有滑动平均部分，上式变为：
$$x_i={\phi }_1{x}_{i-1}+{\phi }_2{x}_{i-2}+...+{\phi }_n{x}_{i-n}+a_i$$
称该表达式为$n$阶自回归模型$AR(n)$。

1.2 模型的参数估计

AR模型的参数估计就是根据已知的观测数据 ，按照某一方法估计出常系数 ${\phi }_i$ 和 ${\sigma }_a^2$ 这 $n+1$ 个参数。参数估计法分为两大类：直接法和间接法。直接法包括最小二乘估计法、解Yuler-Walker方程法等；间接法包括LUD法、Burg法等。由于使用最小二乘法进行参数估计时有着精度高、简单且能够做到无偏估计等优点，因此采用最小二乘法求解AR模型参数。

1.3 模型的定阶方法

AR模型阶次的确定对模型预测的准确性有很大影响，当模型阶数小于最佳模型阶次时，时间序列会产生过于平滑的效果，导致其携带的信息量大大减少；当模型阶数大于最佳模型阶次时，部分噪声依然存在，使得模型拟合的误差增大，影响预测的准确性。
模型定阶通常是根据模型拟合时对原始数据的接近程度来确定利用最佳准则函数，建模时按照准则函数的取值确定模型的优劣。

（1）AIC

赤池信息量准则（Akaike Information Criterion，AIC）是日本学者赤池弘次于1973年提出，并成功应用于AR模型的定阶与选择中。该准则是在熵的概念上而建立的，可以衡量模型拟合的优良性。赤池弘次指出一个时间序列能够被分成若干个局部稳定的时段，每一时段都可用 AR模型来拟合。

AIC准则函数定义为：
$$AIC(k)=2k-2ln(L)$$
其中，$k$是参数的数量，即模型的阶数；$L$为似然函数。

（2）BIC

当训练模型时，增加参数的数量，也就是增加模型的复杂度，从而增大了似然函数，但是模型也会产生过拟合的现象。为了解决这个问题，Schwarz 于1978年提出了贝叶斯信息量准则（Bayes Information Criterion，BIC），又称贝叶斯信息准则。
BIC准则函数定义为：
$$BIC(k)=kln(n)-2ln(L)$$
其中，$k$是参数的数量，即模型的阶数；$n$为样本个数；$L$为似然函数。
$kln(n)$为惩罚项，在样本数量较多的情况下，可有效防止因模型精度过高而造成模型复杂度过高的问题，避免了出现维数灾难的现象。

二、MATLAB代码

1、调用线性预测函数lpc计算模型的参数

clc;
clear;
close all;

randn('state',0); 
noise = randn(50000,1);  % Normalized white Gaussian noise
x = filter(1,[1 1/2 1/3 1/4],noise);
x = x(45904:50000);
a = lpc(x,3);
est_x = filter([0 -a(2:end)],1,x);  % Estimated signal
e = x - est_x;                      % Prediction error
[acs,lags] = xcorr(e,'coeff');      % ACS of prediction error

%   Compare the predicted signal to the original signal
figure;
plot(1:97,x(4001:4097),1:97,est_x(4001:4097),'--');
title('Original Signal vs. LPC Estimate');
xlabel('Sample Number'); ylabel('Amplitude'); grid on;
legend('Original Signal','LPC Estimate')

%   Look at the autocorrelation of the prediction error.
figure;
plot(lags,acs); 
title('Autocorrelation of the Prediction Error');
xlabel('Lags'); ylabel('Normalized Value'); grid on;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

2、利用 Yule-Walker 法计算模型的参数

clc;
clear;
close all;

% Estimate model order using decay of reflection coefficients.
rng default;
y=filter(1,[1 -0.75 0.5],0.2*randn(1024,1));
% Create AR(2) process
[ar_coeffs,NoiseVariance,reflect_coeffs] = aryule(y,10);

% Fit AR(10) model
stem(reflect_coeffs);
axis([-0.05 10.5 -1 1]);
title('Reflection Coefficients by Lag');
xlabel('Lag');ylabel('Reflection Coefficent');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

3、利用 Burg 法计算模型的参数

clc;
clear;
close all;

%  Estimate input noise variance for AR(4) model.
A = [1 -2.7607 3.8106 -2.6535 0.9238]; 
% Generate noise standard deviations
rng default;
noise_stdz = rand(1,50)+0.5;

% Generate column vectors that have corresponding standard deviation
x = bsxfun(@times,noise_stdz,randn(1024,50));

% filter each column using the AR model.
y = filter(1,A,x);

% Compute the estimated coefficients and deviations for each column
[ar_coeffs,NoiseVariance]=arburg(y,4);

% Display the mean value of each estimated polynomial coefficient
estimatedA = mean(ar_coeffs);

% Compare actual vs. estimated variances
plot(noise_stdz.^2,NoiseVariance,'k*');
xlabel('Input Noise Variance');
ylabel('Estimated Noise Variance');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

三、AR相关论文

（1）基于消噪处理岩石声发射信号到达时间的识别方法：在岩石声发射源时差定位的研究中， 信号到达时间是重要信息。岩石声发射信号复杂， 含有大量脉冲干扰与随机噪声， 到达时间可读性差。针对以上问题， 首先对原始声发射信号进行中值滤波和奇异值分解， 消除部分脉冲干扰与随机噪声; 其次进行小波包分解软阈值消噪， 保留信号的主要成分， 提高信噪比， 增强到达时间的可读性; 最后结合信号与噪声的时间序列模型( AutoＲegressive， AＲ 模型) , 第 1 次计算赤池信息量准则的 K 值( Akaike Information Criterion， AIC( K) 值) ， 获得到达时间窗口， 在该窗口内第 2 次计算 AIC( K) 值， 实现了到达时间的自动识别， 避免了在整个信号序列下计算AR模型的阶数与次数。

（2）基于CEEMDAN和AR模型的岩体声发射信号特征提取方法研究：针对岩石声发射信号因非平稳、复杂性特征而难以有效提取特征参数的问题，提出了一种基于自适应完整集成经验分解算法(CEEMDAN)和AR模型的特征提取方法，以红砂岩为研究对象，用CEEMDAN对其破裂阶段的信号进行分解，利用相关系数法选择有效IMF分量，并进行能量归一化处理，以消除其他因素对模型的影响。通过计算各阶段信号AR模型的AIC值，建立最佳阶次模型AR(5)，并通过对比试验验证了模型的准确性，并将提取的５阶AR模型系数作为特征向量，为后续预测岩石失稳现象提供了依据和方法。

参考文献：

[1] 鲜晓东，袁 双，纪松林. 基于消噪处理岩石声发射信号到达时间的识别方法[J]. 煤炭学报, 2015, 40(S1) : 100-106.
[2] 戴聪聪, 罗小燕, 邵凡, 黄小勇, 徐思翔.基于CEEMDAN和AR模型的岩体声发射信号特征提取方法研究[J]. 化工矿物与加工, 2018, 47(06) : 20-24.