最小二乘法求AR模型
阿新 • • 發佈：2019-02-10
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程式碼如下：
對於平穩時間序列而言，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
其中[row,col]代表 row行col列的矩陣，dot代表矩陣點乘運算。
令X的轉置運算為XT，逆矩陣運算為XI。
根據最小二乘的原則，得到引數的計算公式為：
A=(XT dot X)I dot XTdotY
根據該計算公式容易得到p階AR模型引數與預測值計算程式碼：
知道如何計算引數還不夠，還得為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是觀測值與預測值之間殘差的方差。
本例的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的效果更好。
設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程式碼如下：
序列{1978-2014人口死亡率}自相關係數如圖：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;
對於平穩時間序列而言，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
其中[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的效果更好。