ar 65http://www.Sina.com/uto http://www.Sina.com/e gressive模型(自回归模型)利用同一变量以前的表现情况，预测该变量的当前或未来的表现情况。 这种预测方法只与变量本身有关，与其他变量无关，所以称为自回归
数学定义模型：假设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=pN，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. 256.286.346.36.606.906.826.78 ]
6.86.726.646.546.676.706.646.64
6.496.576.56.516.506.466.456.43
6.416.406.426.506.816.937.067.08
7.11 7.14 7.15 7.16 7.16]
A
将mean(x )、var (x )分别作为系列)的平均值和方差，根据自相关系数ACF判断是否为稳定系列。
示例{x}的ACF计算公式为： ACF=(x[I]-mean ) x ) (x[Ik]-mean ) x )、n*var ) x )、0=kN、0=iN-k
python代码如下所示。
导入编号； 导入匹配； #某个k值的ACFdefauto_relate_coef(data，avg，s2，k ) : ef=0. forIinrange(0，len(data )-k ) :ef=ef ) data[I]-AVG ) *(data[i k]-avg )； ef=ef/Len(data )/s2； 返回ef； #计算k0到N-1的所有ACFdefauto_relate_coefs(sample ) : efs=[]； 数据=[ ]； AVG=numpy.mean(sample； S2=numpy.var(sample； array=sample.reshape(1，-1)； forx inarray.flat : data.append (x； forkinrange(0，len ) data ) :ef=auto_relate_coef ) data，avg，s2，k )； EFS.append(ef； return efs； 序列{1978-2014人口死亡率}自相关系数如下。
在稳态时间序列中，ACF系数随着k值的增加而衰减为0的速度快于非稳态随机序列。 由此可见，序列{1978-2014人口死亡率}平稳。
r
根据上述的ar[p]方程式可以得到预测值{y[p]，y[p 1].y[N]}
y [ P1 ]=a [ p ] * x [1] a [ p-1 ] * x [2] . a [1] * x [ p ]
r> 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 XTdotY
根据该计算公式容易得到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的效果更好。