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

Сначала мы создадим массив ежемесячных подсчётов пассажиров авиакомпаний, измеренных тысячами, за период с января 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 просто обозначение для извлечения данных и предоставления их как входных данных любой функции. The controlchart функция также принимает объекты временных рядов непосредственно.

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