AR(Autoregressive)模型（自回归模型）：用同一变量之前的表现情况来预测该变量现在或未来的表现情况，这种预测方法只与变量自己有关，而与其他变量无关，所以称作是自回归。
数学定义模型:假定AR模型是p阶的，对于一组时间序列有观测值{x[1],x[2],.....x[N]}，计算t时刻x的预测值x[t]，其自回归方程：
x[t]=a[1]*x[t-1]+a[2]*x[t-2]....+a[p]*x[t-p]+u[t],1<=p<N,p<=t<=N
其中{a[1],a[2]...a[p]}是对应的参数序列，u[t]是满足N(0,σ^2)的白噪声。
由该数学模型可以看出，AR(p)模型是一种线性预测，由前面的p个x观测值来预测t时刻的值，其本质类似于插值法，其目的都是为了增加有效数据。
AR模型多用于平稳时间序列的预测与拟合，给定一个时间序列，其建模步骤一般如下：
1.判断时间序列是否平稳，可以采用ACF检验、ADF单位根检验等方法。
2若时间序列平稳，则直接转3;若时间序列非平稳，则可采用差分的方法，将其转换为平稳时间序列，转3。
3.计算AR模型的参数（burg算法，最小二乘法，自相关算法等）与定阶（根据AIC准则，SC准则，FPE准则等）。
4.检验3中确定AR模型的拟合度，主要是检验残差序列是否服从N(0,σ^2)白噪声。
5.利用AR模型进行预测。
下面通过实例来分析建模过程：
现有1978-2014年全国人口的死亡率（数据来源于http://www.stats.gov.cn/tjsj/ndsj/）:
[6.25 6.28 6.34 6.36 6.60 6.90 6.82 6.78
6.86 6.72 6.64 6.54 6.67 6.70 6.64 6.64
6.49 6.57 6.56 6.51 6.50 6.46 6.45 6.43
6.41 6.40 6.42 6.50 6.81 6.93 7.06 7.08
7.11 7.14 7.15 7.16 7.16]
1.判断是否为平稳序列
设mean(x),var(x)分别为序列{x}的平均值和方差，根据自相关系数ACF判断是否为平稳序列：
样本{x}的ACF计算公式是：ACF=∑(x[i]-mean(x))*(x[i+k]-mean(x))/(n*var(x))，0<=k<N,0<=i<N-k
python代码如下：
import numpy;
import math;
#计算某一个k值的ACF
def auto_relate_coef(data,avg,s2,k):
ef=0.;
for i in range(0,len(data)-k):
ef=ef+(data[i]-avg)*(data[i+k]-avg);
ef=ef/len(data)/s2;
return ef;
#计算k从0到N-1所有ACF
def auto_relate_coefs(sample):
efs=[];
data=[];
avg=numpy.mean(sample);
s2=numpy.var(sample);
array=sample.reshape(1,-1);
for x in array.flat:
data.append(x);
for k in range(0,len(data)):
ef=auto_relate_coef(data,avg,s2,k);
efs.append(ef);
return efs;
序列{1978-2014人口死亡率}自相关系数如图：
对于平稳时间序列而言，ACF系数随k值的增加衰减到0的速度比非平稳随机序列更快。基于这点可以看出序列{1978-2014人口死亡率}是平稳的。
2.AR模型参数计算与定阶
由上述的AR(p)方程可得到预测值{y[p],y[p+1]....y[N]}
y[p+1]=a[p]*x[1]+a[p-1]*x[2]....a[1]*x[p]
y[p+2]=a[p]*x[2]+a[p-1]*x[3]....a[1]*x[p+1]
.......
y[N]=a[p]*x[N-p]+a[p-1]*x[N-p+1]......a[1]*x[N-1]
将上述方程组写成矩阵形式有：
Y[N-p,1]=X[N-p,p] dotA[p,1]
其中[row,col]代表 row行col列的矩阵，dot代表矩阵点乘运算。
令X的转置运算为XT，逆矩阵运算为XI。
根据最小二乘的原则，得到参数的计算公式为：
A=(XT dot X)I dot XTdot
Y
根据该计算公式容易得到p阶AR模型参数与预测值计算代码：
def ar_least_square(sample,p):
matrix_x=numpy.zeros((sample.size-p,p));
matrix_x=numpy.matrix(matrix_x);
array=sample.reshape(sample.size);
j=0;
for i in range(0,sample.size-p):
matrix_x[i,0:p]=array[j:j+p];
j=j+1;
matrix_y=numpy.array(array[p:sample.size]);
matrix_y=matrix_y.reshape(sample.size-p,1);
matrix_y=numpy.matrix(matrix_y);
#fi为参数序列
fi=numpy.dot(numpy.dot((numpy.dot(matrix_x.T,matrix_x)).I,matrix_x.T),matrix_y);
matrix_y=numpy.dot(matrix_x,fi);
matrix_y=numpy.row_stack((array[0:p].reshape(p,1),matrix_y));
return fi,matrix_y;
知道如何计算参数还不够，还得为AR模型选择一个最优的p值，也就是定阶。
定阶一般步骤为：
(1)确定p值的上限，一般是序列长度N的比例或是lnN的倍数。
(2)在不超过max(p)值的前提下，从1开始根据某一原则确定最优p；
本例中我将p值的上限设为N/2=18，定阶准则用AIC(最小信息准则)和SC(施瓦茨准则)，根据两个准则求得的估计量越小说明阶数越优。
AIC=2*p+N*ln(σ^2) SC=p*ln(N)+N*ln(σ^2)
σ^2是观测值与预测值之间残差的方差。
def ar_aic(rss,p): n=rss.size; s2=numpy.var(rss); return 2*p+n*math.log(s2); def ar_sc(rss,p): n=rss.size; s2=numpy.var(rss); return p*math.log(n)+n*math.log(s2);
本例的AIC和SC：
可以看到当p=18时候，AIC和SC的值均最小，p=19的时候，AIC和SC的值变化比较大。
来看下p=10、18、19时候AR(p)模型的拟合效果(红实线为观测值，蓝虚线为预测值)。
p=10:
p=18
p=19
三幅图可以直观看出p=18时候，AR(18)的拟合效果最好，几乎一模一样。AR(10)虽然效果不如AR(18)，但是扰动在可接受范围内，AR(19)简直丧病，偏离太多。
3.拟合度检验
将AR方程变为下式：
u[t]=x[t]-a[1]*x[t-1]-a[2]*x[t-2]-....-a[p]*x[t-p]
若u[t]是服从N(0,σ^2)的白噪声，则可认为AR(p)是可接受的模型。
本例中用AR(18)计算出的残差u[t],平均值为1.06*10^-6,方差为4.2*10^-4
u[t]的自相关系数如图：
从该图可以看出残差近似服从N(0,σ^2)，因此AR(18)是可以用来拟合和预测的。
总结：
本例采用最小二乘法计算AR模型参数，求得的AR(18)模型效果不俗，缺点在于最小二乘法涉及大量矩阵点乘运算，比较耗时。不止AR模型，还有MA,ARMA,ARIMA模型可以用来拟合和预测平稳时间序列，建模步骤基本一致，相比与AR和MA，ARMA和ARIMA的效果更好。