ar模型的正则方程（matlab中AR模型的计算）

本文目录

• matlab中AR模型的计算
• 诚心求教Matlab关于用Yule-Walker求解AR模型参数问题
• AR 模型的Yule-Walker方程
• ar模型预测误差的方差是什么
• ar模型的传递形式推导
• AR模型的稳定性

matlab中AR模型的计算

=ar(Y,N,Approavh,Window);
%Y观察值
%N AR模型的阶数 取 实整数
%Appproach 表示计算模型的方法 ’fb’,’ls’,’yw’,’burg’,’gl’

诚心求教Matlab关于用Yule-Walker求解AR模型参数问题

function ARMODEL()
Fs = 1000;
t = 0:1/Fs:15;
N = size(t,2) %数据样值点数
randn(’state’,0);
x = cos(2*pi*t*200) + randn(1,N); % 200Hz cosine plus noise
%计算N个取样数据的取样数据自相关函数
rxx = zeros(1,N); %保存取样数据自相关函数的变量
for m = 0:N-1
sum = 0;
for n= 1:N-m
temp1 = x(n)*x(m+n);
sum = sum + temp1;
end
rxx(m+1) = sum/N;
end
%采用Levison-Durbin算法求解AR模型的Yule-Walker模型
%需要确定AR模型理论公式中的参数：白噪声w(n)的方差、方程系数a1……ap（这里包括了模型的阶次）
PMAX = 100; %设定AR模型最高阶次
atemp1 = zeros(1,PMAX+1); %保存方程系数的中间变量
atemp2 = zeros(1,PMAX+1); %保存方程系数的中间变量
deviationtemp1 = zeros; %保存白噪声w(n)方差的中间变量
deviationtemp2 = zeros; %保存白噪声w(n)方差的中间变量
%AR(1)模型：x(n) + a1*x(n-1) = w(n)
%其Yule-Walker方程： R(0)*1 + R(1)*a1 = deviation1;
% R(1)*1 + R(0)*a1 = 0;
%求解方程确定a1、deviation1
atemp1(1) = 1;
atemp1(2) = -rxx(2)/rxx(1);
atemp2 = atemp1;
deviationtemp1 = ( rxx(1)*rxx(1) - rxx(2)*rxx(2) )/rxx(1);
deviationtemp2 = deviationtemp1;
%利用Levison-Durbin迭代算法计算AR模型参数
%根据FPE准则、AIC准则和BIC准则确定AR模型的阶次
%atemp1、deviation1保存第k次的运算结果
%atemp2、deviation2保存第k+1次的运算结果
FPE(1) = deviationtemp1*(N+2)/N;
AIC(1) = log(deviationtemp1) + 2/N;
BIC(1) = log(deviationtemp1) + log(N)/N;
veriance(1) = deviationtemp1;
criteria = 3
for P = 2:PMAX
sum1 = 0;
sum2 = 0;
for i = 2:(P+1)
sum1 = atemp1(i)*rxx(i) + sum1;
end

for i = 1:(P+1)
sum2 = atemp1(i)*rxx(P+2-i) + sum2;
end

deviationtemp1 = rxx(1) + sum1;
dk = sum2;
ref(P) = dk/deviationtemp1;
deviationtemp2 = ( 1 - ref(P)*ref(P) )*deviationtemp1;

for i = 2:(P+1)
atemp2(i) = atemp1(i) - ref(P)*atemp1(P+2-i);
end
%计算AR(P)模型参数
atemp1 = atemp2;
veriance(P) = deviationtemp2

AR 模型的Yule-Walker方程

以AR模型为基础的谱估计由式（1-134）或式（1-135）来计算，这就需要知道模型的阶数p和p个AR系数a_k，以及模型激励源的方差 。为此，必须把这些参数和已知（或估计得到）的自相关函数联系起来，这就是著名的Yule-Walker方程。

Yule-Walker方程可以用两种方法推导。第一种是通过对式（1-134）的P_xx（z）求Z的反变换来得到，另一种是直接由模型的差分方程推导出来，下面介绍第二种推导方法。

将AR模型的差分方程式（1-132）代入x（n）的自相关函数表示式，得

地球物理信息处理基础

地球物理信息处理基础

设AR模型的单位冲激响应是h（n），由方差为 的白噪声序列ε（n）激励h（n）产生输出x（n），于是对于线性时不变系统来说，x（n）=h（n）*ε（n），即

地球物理信息处理基础

如果h（n）是因果的，即m＞0时，h（-m）=0（当 m＜0时，h（m）=0），则上式可写成

地球物理信息处理基础

根据Z变换中的初值定理

地球物理信息处理基础

故上式化为

地球物理信息处理基础

将上式代入式（4-19）

地球物理信息处理基础

得

地球物理信息处理基础

该式称为AR模型的Yule-Walker方程，可以进一步表示为

地球物理信息处理基础

为求取AR模型参数，应先从式（4-20）中选择m＞0的p个方程解出｛a₁，a₂，…，a_p｝，然后代入第一方程（对应于m=0）求出 。自相关函数的头p+1个值是｛r（0），r（1），…，r（p）｝，因此，式（4-20）可表示成下列矩阵形式（这里利用了自相关函数的偶对称性质）：

地球物理信息处理基础

或简单表示为

地球物理信息处理基础

式中：R为x（n）的自相关矩阵；a=［1 a₁a₂… ap］（p个元素）。因此，只要已知或估计出p+1个自相关函数值，即可由该方程解出p+1个模型参数｛a₁，a₂，…，a_p， ｝。

ar模型预测误差的方差是什么

ar(2)模型的方差公式为:Rj=a1R(j-1)+a2R(j-2)。在用 AR 模型 对数据进行建模时AR模型的最佳预测 所谓的最佳预测,即是预测误差的方差为最小意义下的预测预测误差方差: 注意:因为在预测原点对时刻的序列值做出预测时,时刻的序列值并未发生,所以需要首先计算作为时刻序列值的近似替换。
2. 预测误差 2步预测误差方差AR(1)过程模拟的R实现 X_n=0.5X_{n-1}+\epsilon_n, \epsilon_n\sim N(0,1) wn《-
2. X_n=1.5X_{n-1}-0.3X_{n-2}+\epsilon_n 的模拟 》 y《-arima.sim(model=list(ar=c(
3. 模拟模型 X_n=0.5X_{n-1}+0.3X_{

ar模型的传递形式推导

2.AR(p)过程

\{X_n:n=0,\pm1,\pm2,\cdots\} 为一时间序列，满足对 \forall n \in ZX_n=a_0+a_1X_{n-1}+a_2X_{n-2}+\cdots+a_pX_{n-p}+\epsilon_n

\{\epsilon_n:n=0,\pm1,\pm2,\cdots\} 为白噪声序列， Var(\epsilon_n)=\sigma^2 。称 \{X_n\} 为 p 阶自回归过程，记为AR(p)。

3.特征方程

1-a_1x-a_2x^2-\cdots-a_px^p=0

结论：AR(p)过程是平稳的当且仅当特征方程解的模大于1.

4.均值为0的AR(p)过程

令 \mu 为平稳的AR(p)过程均值，那么 X_n-\mu=a_1(X_{n-1}-\mu)+a_2(X_{n-2}-\mu)+\cdots+a_k(X_{n-p}-\mu)+\epsilon_n

因此在以下的讨论中，我们假定AR(p)过程为

X_n=a_1X_{n-1}+a_2X_{n-2}+\cdots+a_kX_{n-p}+\epsilon_n

5.Yule－Walker方程组

令 \rho_i,i=1,2,\cdots,p 为 i 阶自相关系数，则

\left\{ \begin{aligned} \rho_1 &= a_1 +a_2\rho_2+a_3\rho_3+\cdots+a_p\rho_p\\ \rho_2 & =a_1 \rho_1+a_2+a_3\rho_3+\cdots+a_p\rho_p\\ &\cdots&\\ \rho_p & =a_1 \rho_1+a_2\rho_2+a_3\rho_3+\cdots+a_p \end{aligned} \right.

解方程组可得自相关系数 \rho_i,i=1,2,\cdots,p 。

然后利用递推关系求出各阶自相关系数：

\rho_k=a_1 \rho_{k-1}+a_2\rho_{k-2}+a_3\rho_{k-3}+\cdots+a_p\rho_{k-p},k\geq p+1 .

6. 平稳的AR(1) 过程

X_n=a_1X_{n-1}+\epsilon_n , a_1\ne 0

性质

(1) |a_1|《1 ;

(2) Var(X_n)=\frac{\sigma^2}{1-a^2_1} ；

(3) 自相关系数 \rho_k=a^k_1,k\geq0 .

AR模型的稳定性

AR（p）模型稳定的充分必要条件是H（z）的极点（即A（z）的根）都在单位圆内。如果Yule-Walker方程的系数矩阵是正定的，则其解a_k（k=1，2，…，p）所构成的A（z）的根都在单位圆内。在用Levinson算法进行递推计算的过程中，还可得到各阶AR模型激励信号的方差 （k=1，2，…，p），它们都应当是大于零的，即 。根据式（4-25）可知，必有｜γ_k+1｜＜1和 （k=1，2，…，p）。这就是说，在Levinson算法递推计算过程中，如果有 ＜ 或，｜γ_k+1｜＜1，则AR（p）模型一定是稳定的。反之，稳定的AR（p）模型将具有以下性质：

（1）H（z）的全部极点或A（z）的所有根都在单位圆内。

（2）自相关矩阵是正定的。

（3）激励信号的方差（能量）随阶次增加而递减，即 。

（4）反射系数的模恒小于1，即，γ_k＜1，k=1，2，…，p。

但在实际应用中，Levinson算法的已知数据（自相关值）是由x_N（n）来估计的，有限字长效应有可能造成大的误差，致使估计出来的AR（p）参数所构成的A（z）的根跑到单位圆上或外，从而使模型失去稳定。在递推计算过程中如果出现这种情况，将导致 或，｜γ_k｜≥1，即停止递推计算。

若将式（4-22）中的自相关矩阵R定为

地球物理信息处理基础

并记其行列式的值为detR_p+1。矩阵R_p+1与AR（p）模型稳定性的关系有以下三个结论。

结论1：如果R_p+1是正定的，那么，由Yule-Walker方程解出的a_k（k=1，2，…，p）构成的p阶AR模型是稳定的，且是唯一的。也即A（z）的零点都在单位圆内。此性质为AR模型的最小相位性质。

结论2：若x（n）由p个复正弦波组成，即

地球物理信息处理基础

式中：C_k、ω_k为常数；φ_k是在（-π，π）内均匀分布的零均值随机变量；x（n）的自相关函数为

地球物理信息处理基础

则由前p+1个值r_xx（0）、r_xx（1）、…、r_xx（p）组成的自相关矩阵R_p+1是奇异的，而R₁、R₂、…、R_p是正定的，即

det R_p+1=0，det R_k＞0 k=1，2，…，p （4-55）

结论2说明，一般情况下，若x（n）由复正弦波组成，R_M是其M×M的自相关矩阵，那么，当M＞p时，R_M的秩最大为p，即rankR_M=p，但若x（n）由p个实正弦波组成，则R_M的秩最大为2p。

结论3：若x（n）由p个正弦波组成（实的或复的），则x（n）是完全可以预测的，即预测误差为零。

结论2给出了R_p+1何时奇异、何时正定的条件，它和结论3一起揭示了正弦信号的某些性质。特别需要指出的是：用AR模型对纯正弦信号建模是不合适的，可能会出现自相关矩阵为奇异的情况。但是，在信号处理中经常要用正弦信号作为实验信号以检验某个算法或系统的性能。为了克服自相关矩阵奇异的情况，最常用的方法是给正弦信号加上白噪声，这样det R_p+1不会等于零。