A moving average with different functions and varying time-frames

Question

I have a matrix time-series data for 8 variables with about 2500 points (~10 years of mon-fri) and would like to calculate the mean, variance, skewness and kurtosis on a 'moving average' basis.

Lets say frames = [100 252 504 756] - I would like calculate the four functions above on over each of the (time-)frames, on a daily basis - so the return for day 300 in the case with 100 day-frame, would be [mean variance skewness kurtosis] from the period day201-day300 (100 days in total)... and so on.

I know this means I would get an array output, and the the first frame number of days would be NaNs, but I can't figure out the required indexing to get this done...

Answer 1

This is an interesting question because I think the optimal solution is different for the mean than it is for the other sample statistics.

I've provided a simulation example below that you can work through.

First, choose some arbitrary parameters and simulate some data:

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

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

For the mean, use filter to obtain a moving average:

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

I had originally thought to solve this problem using conv as follows:

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

But as @PhilGoddard pointed out in the comments, the filter approach avoids the need for the loop.

Also note that I've chosen to make the dates in the output matrix correspond to the dates in X so in later work you can use the same subscripts for both. Thus, the first WindowLength-1 observations in MeanMA will be nan.

For the variance, I can't see how to use either filter or conv or even a running sum to make things more efficient, so instead I perform the calculation manually at each iteration:

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

We could speed things up slightly by exploiting the fact that we have already calculated the mean moving average. Simply replace the within loop line in the above with:

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

However, I doubt this will make much difference.

If anyone else can see a clever way to use filter or conv to get the moving window variance I'd be very interested to see it.

I leave the case of skewness and kurtosis to the OP, since they are essentially just the same as the variance example, but with the appropriate function.

A final point: if you were converting the above into a general function, you could pass in an anonymous function as one of the arguments, then you would have a moving average routine that works for arbitrary choice of transformations.

Final, final point: For a sequence of window lengths, simply loop over the entire code block for each window length.

Answer 2

I have managed to produce a solution, which only uses basic functions within MATLAB and can also be expanded to include other functions, (for finance: e.g. a moving Sharpe Ratio, or a moving Sortino Ratio). The code below shows this and contains hopefully sufficient commentary.

I am using a time series of Hedge Fund data, with ca. 10 years worth of daily returns (which were checked to be stationary - not shown in the code). Unfortunately I haven't got the corresponding dates in the example so the x-axis in the plots would be 'no. of days'.

% 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')