В этом примере показано, как визуализировать и анализировать данные временных рядов с помощью timeseries возразите и regress функция.
Сначала мы создаем массив ежемесячных количеств авиапассажиров, измеренных в тысячах, в течение января 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™.