我正在尝试创建一个程序,使用Hannan-Rissanen算法计算ARMA(p,q)自回归移动平均随机过程的样本参数。
我遇到困难的主要步骤是计算时间序列的自协方差函数。
该程序应该采用n×1维列向量Y并计算由下式给出的k×1维列向量γ^ hat:
其中Ybar是Y元素的平均值。
如何有效地计算上述总和? (显然for循环可以工作,但我正试图在矢量化numpy操作上做得更好)因为我正在使用它作为一种学习经验,我宁愿不使用除了像np.sum
或np.mean
这样的非常基本的函数之外的任何numpy函数。
以下类似的问题已被提出,但不能完全回答我的问题:
Computing autocorrelation of vectors with numpy(使用np.correlate
)
(其他一些人遇到了使用更高级的numpy
函数的相同问题,或者没有像我希望的那样吐出向量。)
这是替换np.correlate
的一种方法(我认为这是主要的困难;我也假设你不想手工编写fft):
def autocorr_direct(a, debug=False):
n, _ = a.shape
out = np.zeros((n+1, 2*n-1), a.dtype)
out.reshape(-1)[:2*n*n].reshape(n, 2*n)[::-1, :n] = a*a.T
if debug:
print(out.reshape(-1)[:2*n*n].reshape(n, 2*n))
print(out)
return out.sum(0)
例如:
>>> a = np.array([[1, 1, 2, -1]]).T
>>> autocorr_direct(a, True)
[[-1 -1 -2 1 0 0 0 0]
[ 2 2 4 -2 0 0 0 0]
[ 1 1 2 -1 0 0 0 0]
[ 1 1 2 -1 0 0 0 0]]
[[-1 -1 -2 1 0 0 0]
[ 0 2 2 4 -2 0 0]
[ 0 0 1 1 2 -1 0]
[ 0 0 0 1 1 2 -1]
[ 0 0 0 0 0 0 0]]
array([-1, 1, 1, 7, 1, 1, -1])
>>> np.correlate(a[:, 0], a[:, 0], 'full')
array([-1, 1, 1, 7, 1, 1, -1])
注意剪切方形阵列a[::-1]*a.T
的重塑技巧。
笔记2;从1D矢量X
获得列向量使用X[:, None]
。