Линейная регрессия в частотной области

В этом примере показано, как использовать дискретное преобразование Фурье для построения модели линейной регрессии для временного ряда. Временные ряды, используемые в этом примере, - это месячное число случайных смертей в США с 1973 по 1979 год. Данные опубликованы в Brockwell and Davis (2006). Первоисточник - Национальный совет по безопасности США.

Введите данные. Копировать exdata в рабочую область MATLAB ®.

exdata = [
        9007        7750        8162        7717        7792        7836
        8106        6981        7306        7461        6957        6892
        8928        8038        8124        7776        7726        7791
        9137        8422        7870        7925        8106        8129
       10017        8714        9387        8634        8890        9115
       10826        9512        9556        8945        9299        9434
       11317       10120       10093       10078       10625       10484
       10744        9823        9620        9179        9302        9827
        9713        8743        8285        8037        8314        9110
        9938        9129        8433        8488        8850        9070
        9161        8710        8160        7874        8265        8633
        8927        8680        8034        8647        8796        9240];

exdata является матрицей 12 на 6. Каждый столбец exdata содержит данные за 12 месяцев. Первая строка каждого столбца содержит количество случайных смертей в США за январь соответствующего года. Последняя строка каждого столбца содержит количество случайных смертей в США за декабрь соответствующего года.

Преобразуйте матрицу данных в временной ряд 72 на 1 и постройте график данных за 1973-1978 годы.

ts = reshape(exdata,72,1);
years = linspace(1973,1979,72);

plot(years,ts,'o-','MarkerFaceColor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')

Визуальный осмотр данных показывает, что количество случайных смертей варьируется периодически. Период колебания, по-видимому, составляет примерно 1 год (12 месяцев). Периодический характер данных позволяет предположить, что соответствующая модель может быть

где $$\lambda$$ - общее среднее значение, $$N$$ - длина временного ряда, и $$\lambda (n$$) - последовательность белого шума независимых и идентично распределенных (IID) гауссовых случайных величин с нулевым средним и некоторой дисперсией. Термин аддитивный шум учитывает случайность, присущую данным. Параметрами модели являются общее среднее и амплитуды косинусов и синусов. Модель линейна в параметрах.

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

Выполните спектральный анализ данных, чтобы определить, какие частоты вносят значительный вклад в изменчивость данных. Поскольку общее среднее значение сигнала составляет приблизительно 9000 и пропорционально преобразованию Фурье на частоте 0, вычитайте среднее значение до спектрального анализа. Это уменьшает коэффициент Фурье большой величины при частоте 0 и облегчает обнаружение любых значительных колебаний. Частоты в преобразовании Фурье разнесены с интервалом, который является обратным от длины временного ряда, 1/72. Выборка данных ежемесячно, наибольшая частота в спектральном анализе - 1 цикл/2 месяца. В этом случае удобно смотреть на спектральный анализ в терминах циклов/год, так что масштабировать частоты соответственно для визуализации.

tsdft = fft(ts-mean(ts));
freq = 0:1/72:1/2;

plot(freq.*12,abs(tsdft(1:length(ts)/2+1)),'o-', ...
    'MarkerFaceColor','auto')
xlabel('Cycles/Year')
ylabel('Magnitude')
ax = gca;
ax.XTick = [1/6 1 2 3 4 5 6];

Исходя из величин, частота 1 цикл/12 месяцев является наиболее значительным колебанием в данных. Величина за 1 цикл/12 месяцев более чем в два раза превышает любую другую величину. Однако спектральный анализ показывает, что в данных имеются и другие периодические компоненты. Например, по-видимому, существуют периодические компоненты на гармониках (целых кратных) 1 цикл/12 месяцев. Как представляется, существует также периодический компонент с периодом в 1 цикл/72 месяца.

На основе спектрального анализа данных подгоняют простую модель линейной регрессии с использованием косинусного и синусоидального члена с частотой наиболее значимой составляющей: 1 цикл/год (1 цикл/12 месяцев).

Определите частотный блок в дискретном преобразовании Фурье, который соответствует 1 циклу/12 месяцам. Поскольку частоты разнесены на 1/72, а первый бункер соответствует 0 частоте, правильный бункер равен 72/12 + 1. Это частотный блок положительной частоты. Необходимо также включить ячейку частоты, соответствующую отрицательной частоте: -1 цикл/12 месяцев. При индексировании MATLAB частотный блок отрицательной частоты равен 72-72/12 + 1.

Создайте вектор нулей 72 на 1. Заполните соответствующие элементы вектора коэффициентами Фурье, соответствующими положительной и отрицательной частоте 1 цикл/12 месяцев. Инвертировать преобразование Фурье и добавить общее среднее значение, чтобы получить соответствие с данными о случайной смерти.

freqbin = 72/12;
freqbins = [freqbin 72-freqbin]+1;
tsfit = zeros(72,1);
tsfit(freqbins) = tsdft(freqbins);
tsfit = ifft(tsfit);
mu = mean(ts);
tsfit = mu+tsfit;

Постройте график исходных данных вместе с подобранными сериями, используя два коэффициента Фурье.

plot(years,ts,'o-','MarkerFaceColor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')
hold on
plot(years,tsfit,'linewidth',2)
legend('Data','Fitted Model')
hold off

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

Чтобы оценить, насколько адекватно одна частота 1 цикл/12 месяцев учитывает наблюдаемые временные ряды, сформируйте остатки. Если остатки напоминают последовательность белого шума, простая линейная модель с одной частотой адекватно смоделировала временной ряд.

Чтобы оценить остатки, используйте автокорреляционную последовательность с 95% -ными интервалами достоверности для белого шума.

resid = ts-tsfit;
[xc,lags] = xcorr(resid,50,'coeff');

stem(lags(51:end),xc(51:end),'filled')
hold on
lconf = -1.96*ones(51,1)/sqrt(72);
uconf = 1.96*ones(51,1)/sqrt(72);
plot(lags(51:end),lconf,'r')
plot(lags(51:end),uconf,'r')
xlabel('Lag')
ylabel('Correlation Coefficient')
title('Autocorrelation of Residuals')
hold off

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

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

tsfit2dft = zeros(72,1);
[Y,I] = sort(abs(tsdft),'descend');
indices = I(1:6);
tsfit2dft(indices) = tsdft(indices);

Продемонстрируйте, что сохранение только 6 из 72 коэффициентов Фурье (3 частоты) сохраняет большую часть энергии сигнала. Во-первых, продемонстрировать, что сохранение всех коэффициентов Фурье дает эквивалентность энергии между исходным сигналом и преобразованием Фурье.

norm(1/sqrt(72)*tsdft,2)/norm(ts-mean(ts),2)

ans = 1.0000

Отношение равно 1. Теперь изучите отношение энергии, где сохраняются только 3 частоты.

norm(1/sqrt(72)*tsfit2dft,2)/norm(ts-mean(ts),2)

ans = 0.8991

Почти 90% энергии сохраняется. Эквивалентно, 90% дисперсии временного ряда учитывается 3 частотными компонентами.

Сформировать оценку данных на основе 3 частотных компонентов. Сравните исходные данные, модель с одной частотой и модель с 3 частотами.

tsfit2 = mu+ifft(tsfit2dft,'symmetric');

plot(years,ts,'o-','markerfacecolor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')
hold on
plot(years,tsfit,'linewidth',2)
plot(years,tsfit2,'linewidth',2)
legend('Data','1 Frequency','3 Frequencies')
hold off

Использование 3 частот улучшило соответствие исходному сигналу. Это можно увидеть, изучив автокорреляцию остатков из 3-частотной модели.

resid = ts-tsfit2;
[xc,lags] = xcorr(resid,50,'coeff');

stem(lags(51:end),xc(51:end),'filled')
hold on
lconf = -1.96*ones(51,1)/sqrt(72);
uconf = 1.96*ones(51,1)/sqrt(72);
plot(lags(51:end),lconf,'r')
plot(lags(51:end),uconf,'r')
xlabel('Lag')
ylabel('Correlation Coefficient')
title('Autocorrelation of Residuals')
hold off

Использование 3 частот привело к образованию остатков, которые более близко приближаются к процессу белого шума.

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

X = ones(72,7);
X(:,2) = cos(2*pi/72*(0:71))';
X(:,3) = sin(2*pi/72*(0:71))';
X(:,4) = cos(2*pi*6/72*(0:71))';
X(:,5) = sin(2*pi*6/72*(0:71))';
X(:,6) = cos(2*pi*12/72*(0:71))';
X(:,7) = sin(2*pi*12/72*(0:71))';
beta = X\ts;
tsfit_lm = X*beta;
max(abs(tsfit_lm-tsfit2))

ans = 1.2733e-11

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

Для получения более подробной информации о спектральном анализе во временных рядах и эквивалентности регрессии во временной области см. (Shumway and Stoffer, 2006).

Хотя спектральный анализ может ответить, какие периодические компоненты вносят значительный вклад в изменчивость данных, он не объясняет, почему эти компоненты присутствуют. При тщательном изучении этих данных можно увидеть, что минимальные значения в 12-месячном цикле, как правило, приходятся на февраль, а максимальные - на июль. Правдоподобным объяснением этих данных является то, что люди, естественно, более активны летом, чем зимой. К сожалению, в результате этой повышенной активности повышается вероятность возникновения несчастных случаев со смертельным исходом.

Ссылки

Броквелл, Питер Дж. и Ричард А. Дэвис. Временные ряды: теория и методы. Нью-Йорк: Спрингер, 2006.

Шумвей, Роберт Х. и Дэвид С. Стоффер. Анализ временных рядов и его применения с помощью примеров R. Нью-Йорк: Спрингер, 2006.

См. также

fft | ifft | xcorr