不同函数和不同时间框架的移动平均线

Question

我有一个包含8个变量的矩阵时间序列数据，大约有2500个数据点（约为10年的周一至周五），我想按“移动平均”的方式计算平均数、方差、偏度和峰度。例如，假设frame = [100 252 504 756]，我想针对每个（时间）帧上的日常基础进行以上四个函数的计算 - 因此，在具有100天帧的情况下，第300天的返回值将是从day201-day300（总共100天）期间的[mean variance skewness kurtosis]......以此类推。

我知道这意味着我会得到一个数组输出，并且前frame天的数值将是NaN，但是我无法找出所需的索引来完成此操作...

- n1k31t4

Answer 1

这是一个有趣的问题，因为我认为对于平均值来说最优解与其他样本统计量不同。下面提供了一个模拟示例供您参考。

首先，选择一些任意参数并模拟一些数据:

%#Set some arbitrary parameters
T = 100; N = 5;
WindowLength = 10;

%#Simulate some data
X = randn(T, N);

对于均值，使用filter来获取移动平均值：

MeanMA = filter(ones(1, WindowLength) / WindowLength, 1, X);
MeanMA(1:WindowLength-1, :) = nan;

我最初想使用 conv 来解决这个问题，如下所示：

MeanMA = nan(T, N);
for n = 1:N
    MeanMA(WindowLength:T, n) = conv(X(:, n), ones(WindowLength, 1), 'valid');
end
MeanMA = (1/WindowLength) * MeanMA;

但是，正如@PhilGoddard在评论中指出的那样，filter方法避免了使用循环的必要。

另外请注意，我选择使输出矩阵中的日期与X中的日期对应，因此，在后续工作中您可以为两者使用相同的下标。 因此，MeanMA中前WindowLength-1个观测值将为nan。

至于方差，我无法想到如何使用filter或conv甚至是一个运行总和来实现更高效的计算，因此我在每次迭代时手动执行计算：

VarianceMA = nan(T, N);
for t = WindowLength:T
    VarianceMA(t, :) = var(X(t-WindowLength+1:t, :));
end

我们可以利用已经计算出来的平均移动平均值，稍微加快速度。只需将上述代码中的循环内部行替换为：

VarianceMA(t, :) = (1/(WindowLength-1)) * sum((bsxfun(@minus, X(t-WindowLength+1:t, :), MeanMA(t, :))).^2);

不过，我怀疑这样做并不会有太大的区别。

如果有人能想到一个巧妙地使用 filter 或者 conv 来获得移动窗口方差的方法，我将非常感兴趣。

至于偏度和峰度的情况，我将其留给OP，因为它们与方差示例本质上是相同的，只需适当的函数即可。

最后要说的一点是：如果您将上述内容转换为通用函数，可以将匿名函数作为参数之一传入，那么您将拥有一个适用于任意变换选择的移动平均计算程序。

最后一个观点：对于一系列的窗口长度，只需为每个窗口长度循环整个代码块即可。

- Colin T Bowers

---

对于平均值，最好使用filter函数，这样可以消除循环的需要。对于其他统计数据，只能使用循环。 - Phil Goddard

@PhilGoddard 谢谢。直到现在，filter函数在我的视线范围之外。我已经更新了我的答案，包括这种方法。 - Colin T Bowers

谢谢您的建议 - 它很有帮助。然而，我不认为它必须如此复杂 - 请看下面我的解决方案。也许您的解决方案以后可以做更多的事情？？我的解决方案至少解决了我的初始问题，并且可以扩展以包括更多功能。 - n1k31t4

@DexterMorgan 你知道吗，你的均值/方差/偏度/峰度解决方案与我的方差解决方案完全相同吗？所以当你说“复杂”时，我认为你是指我的平均值解决方案——具体来说，我的使用filter函数。对于Matlab的一个通用规则（尽管这些天由于JIT编译器的进步而不那么重要）是通过向量化代码来避免循环。这就是我答案中filter函数的目的。 - Colin T Bowers

是的，对于平均值来说，过滤函数确实更好 - 但我想对几个不同的函数进行操作，而不仅仅是平均值。我发表了我的答案，因为它对我有用，我认为它也可能帮助其他人。 - n1k31t4

@DexterMorgan 没问题，但我恐怕不能给它点赞，因为它本质上与我在答案中提到的方差是一样的。 - Colin T Bowers

Answer 2

我已经成功地开发了一个解决方案，它只使用了MATLAB中的基本函数，并且还可以扩展以包括其他函数（例如金融：移动夏普比率或移动索提诺比率）。下面的代码显示了这个解决方案，并包含了足够的注释。

我正在使用一个对冲基金数据的时间序列，大约有10年的每日回报率（已检查为平稳-未在代码中显示）。不幸的是，在这个例子中我没有相应的日期，所以图表中的x轴将是“天数”。

% start by importing the data you need - here it is a selection out of an
% excel spreadsheet
returnsHF =  xlsread('HFRXIndices_Final.xlsx','EquityHedgeMarketNeutral','D1:D2742');

% two years to be used for the moving average. (250 business days in one year)
window = 500;

% create zero-matrices to fill with the MA values at each point in time. 
mean_avg = zeros(length(returnsHF)-window,1);
st_dev = zeros(length(returnsHF)-window,1);
skew = zeros(length(returnsHF)-window,1);
kurt = zeros(length(returnsHF)-window,1);

% Now work through the time-series with each of the functions (one can add
% any other functions required), assinging the values to the zero-matrices
for count = window:length(returnsHF)

% This is the most tricky part of the script, the indexing in this section
% The TwoYearReturn is what is shifted along one period at a time with the
% for-loop. 
TwoYearReturn = returnsHF(count-window+1:count);
mean_avg(count-window+1) = mean(TwoYearReturn);
st_dev(count-window+1) = std(TwoYearReturn);
skew(count-window+1) = skewness(TwoYearReturn);
kurt(count-window +1) = kurtosis(TwoYearReturn);
end

% Plot the MAs
subplot(4,1,1), plot(mean_avg)
title('2yr mean')
subplot(4,1,2), plot(st_dev)
title('2yr stdv')
subplot(4,1,3), plot(skew)
title('2yr skewness')
subplot(4,1,4), plot(kurt)
title('2yr kurtosis')

- n1k31t4