最小二乘法求AR模型

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;

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;

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);

可以看到当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的效果更好。