Данные о воздушных пассажирах

Сначала мы создаем массив ежемесячных подсчетов пассажиров авиакомпании, измеренных тысячами, за период с января 1949 года по декабрь 1960 года.

%   1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960
y = [112  115  145  171  196  204  242  284  315  340  360  417    % Jan
     118  126  150  180  196  188  233  277  301  318  342  391    % Feb
     132  141  178  193  236  235  267  317  356  362  406  419    % Mar
     129  135  163  181  235  227  269  313  348  348  396  461    % Apr
     121  125  172  183  229  234  270  318  355  363  420  472    % May
     135  149  178  218  243  264  315  374  422  435  472  535    % Jun
     148  170  199  230  264  302  364  413  465  491  548  622    % Jul
     148  170  199  242  272  293  347  405  467  505  559  606    % Aug
     136  158  184  209  237  259  312  355  404  404  463  508    % Sep
     119  133  162  191  211  229  274  306  347  359  407  461    % Oct
     104  114  146  172  180  203  237  271  305  310  362  390    % Nov
     118  140  166  194  201  229  278  306  336  337  405  432 ]; % Dec
% Source:
% Hyndman, R.J., Time Series Data Library,
% http://www-personal.buseco.monash.edu.au/~hyndman/TSDL/.
% Copied in October, 2005.

Создать объект временного ряда

При создании объекта временного ряда можно сохранить информацию о времени вместе со значениями данных. У нас есть ежемесячные данные, поэтому мы создаем массив дат и используем его вместе с данными Y для создания объекта временных рядов.

yr = repmat((1949:1960),12,1);
mo = repmat((1:12)',1,12);
time = datestr(datenum(yr(:),mo(:),1));
ts = timeseries(y(:),time,'name','AirlinePassengers');
ts.TimeInfo.Format = 'dd-mmm-yyyy';
tscol = tscollection(ts);
plot(ts)

Изучение тренда и сезонности

Эта серия, по-видимому, имеет сильную сезонную составляющую, с тенденцией, которая может быть линейной или квадратичной. Кроме того, величина сезонных колебаний увеличивается по мере увеличения общего уровня. Возможно, преобразование логарифма сделает сезонные колебания более постоянными. Сначала мы изменим масштаб оси.

h_gca = gca;
h_gca.YScale = 'log';

Представляется, что было бы проще смоделировать сезонный компонент в логарифмической шкале. Мы создадим новый временной ряд с преобразованием журнала.

tscol = addts(tscol,log(ts.data),'logAirlinePassengers');
logts = tscol.logAirlinePassengers;

Теперь нарисуем средние значения за год с наложением ежемесячных отклонений. Мы хотим увидеть, является ли изменение от месяца к месяцу постоянным в течение лет. Для таких манипуляций, рассматривающих данные как матрицу в ежемесячном формате, удобнее работать с исходной матрицей данных.

t = reshape(datenum(time),12,12);
logy = log(y);
ymean = repmat(mean(logy),12,1);
ydiff = logy - ymean;
x = yr + (mo-1)/12;
plot(x,ymean,'b-',x,ymean+ydiff,'r-')
title('Monthly variation within year')
xlabel('Year')

Теперь давайте обратим вспять годы и месяцы, и попробуем посмотреть, является ли тренд год к году постоянным для каждого месяца.

h_gca = gca;
h_gca.Position = [0.13 0.58 0.78 0.34];
subplot(2,1,2);
t = reshape(datenum(time),12,12);
mmean = repmat(mean(logy,2),1,12);
mdiff = logy - mmean;
x = mo + (yr-min(yr(:)))/12;
plot(x',mmean','b-',x',(mmean+mdiff)','r-')
title('Yearly trend within month')
xlabel('Month')

Модельный тренд и сезонность

Попробуем смоделировать эту серию как линейный тренд плюс сезонная составляющая.

subplot(1,1,1);
X = [dummyvar(mo(:)) logts.time];
[b,bint,resid] = regress(logts.data,X);
tscol = addts(tscol,X*b,'Fit1')

Time Series Collection Object: unnamed

Time vector characteristics

      Start date            01-Jan-1949
      End date              01-Dec-1960

Member Time Series Objects:

      AirlinePassengers
      logAirlinePassengers
      Fit1

plot(logts)
hold on
plot(tscol.Fit1,'Color','r')
hold off
legend('Data','Fit','location','NW')

На основе этого графика посадка представляется хорошей. Различия между фактическими данными и соответствующими значениями вполне могут быть достаточно малыми для наших целей.

Но давайте попробуем исследовать это еще. Мы хотели бы, чтобы остатки выглядели независимыми. Если есть автокорреляция (корреляция между соседними остатками), то может быть возможность смоделировать это и сделать нашу подгонку лучше. Давайте создадим временной ряд из остатков и построим его график.

tscol = addts(tscol,resid,'Resid1');
plot(tscol.Resid1)

Остатки не выглядят независимыми. На самом деле корреляция между соседними остатками выглядит довольно сильной. Мы можем проверить это формально, используя тест Дурбина-Уотсона.

[p,dw] = dwtest(tscol.Resid1.data,X)

p = 7.7787e-30

dw = 0.4256

Низкое значение p для статистики Дурбина-Уотсона является показателем того, что остатки коррелируются во времени. Типичным ограничением для тестов гипотез является определение того, что p < 0,05 является значимым. Здесь очень небольшое значение p дает убедительные доказательства того, что остатки коррелированы.

Мы можем попытаться изменить модель, чтобы удалить автокорреляцию. Общая форма кривой высока в середине и низка на концах. Это говорит о том, что мы должны допустить квадратичный тренд. Тем не менее, также кажется, что автокорреляция останется после того, как мы добавим этот термин. Попробуем.

X = [dummyvar(mo(:)) logts.time logts.time.^2];
[b2,bint,resid2] = regress(logts.data,X);
tscol = addts(tscol,resid2,'Resid2');
plot(tscol.Resid2)

[p,dw] = dwtest(tscol.Resid2.data,X)

p = 8.7866e-20

dw = 0.6487

Сложение квадратного члена действительно удалило выраженную кривизну в исходном остаточном графике, но и график, и новый тест Дурбина-Уотсона показывают, что в остатках все еще есть значительная корреляция.

Такая автокорреляция может быть результатом других причин, которые не зафиксированы в нашей переменной X. Возможно, мы могли бы собрать другие данные, которые помогли бы нам улучшить нашу модель и уменьшить корреляцию. В отсутствие других данных мы можем просто добавить другой параметр в модель, чтобы представить автокорреляцию. Давайте сделаем это, удалив квадрат и используя авторегрессионную модель для ошибки.

В авторегрессионном процессе у нас есть две стадии:

   Y(t) = X(t,:)*b + r(t)       % regression model for original data
   r(t) = rho * r(t-1) + u(t)   % autoregressive model for residuals

В отличие от обычной регрессионной модели, когда мы хотели бы остаточный ряд r(t) чтобы быть набором независимых значений, эта модель позволяет остаткам следовать авторегрессивной модели с собственным термином ошибки u(t) который состоит из независимых значений.

Чтобы создать эту модель, мы хотим написать анонимную функцию f для вычисления соответствующих значений Yfit, так что Y-Yfit дает значения u:

   Yfit(t) = rho*Y(t-1) + (X(t,:) - rho*X(t-1,:))*b

В этой анонимной функции мы объединяем [rho; b] в один вектор параметров c. Полученные остатки выглядят гораздо ближе к некоррелированному ряду.

r = corr(resid(1:end-1),resid(2:end));  % initial guess for rho
X = [dummyvar(mo(:)) logts.time];
Y = logts.data;
f = @(c,x) [Y(1); c(1)*Y(1:end-1) + (x(2:end,:)- c(1)*x(1:end-1,:))*c(2:end)];
c = nlinfit(X,Y,f,[r; b]);

u = Y - f(c,X);
tscol = addts(tscol,u,'ResidU');
plot(tscol.ResidU);

Резюме

В этом примере представлена иллюстрация использования объекта временных рядов MATLAB ® вместе с элементами из панели инструментов статистики и машинного обучения. Это просто использовать ts.data запись для извлечения данных и передачи их в качестве входных данных для любой функции. controlchart функция также принимает объекты временных рядов напрямую.

Более детальный анализ возможен с использованием функций, специально разработанных для временных рядов, таких как в Econometrics Toolbox™ и System Identification Toolbox™.