Данные авиапассажира

Сначала мы создаем массив ежемесячных количеств авиапассажиров, измеренных в тысячах, в течение января 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® наряду с функциями из Statistics and Machine Learning Toolbox. Это просто в использовании ts.data обозначение, чтобы извлечь данные и предоставить его, как введено к любой функции. controlchart функция также принимает объекты временных рядов непосредственно.

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