В этом примере показано, как использовать дискретное преобразование Фурье для создания линейной регрессионой модели для временных рядов. Временные ряды, используемые в этом примере, являются ежемесячным количеством случайных смертей в Соединенных Штатах с 1973 по 1979 год. Данные опубликованы в Броквелле и Дэвисе (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 месяцев). Периодический характер данных предполагает, что соответствующая модель может быть
где - общее среднее, - длина временных рядов, и является белой шумовой последовательностью независимых и одинаково распределенных (IID) случайных переменных Гауссова с нулем среднего и некоторым отклонением. Термин аддитивного шума учитывает случайность, присущую данным. Параметрами модели являются общее среднее и амплитуды косинусов и синусов. Модель линейна в параметрах.
Чтобы создать линейную регрессионую модель в временной интервал, вы должны задать, какие частоты использовать для косинусов и синусов, сформировать матрицу проекта и решить нормальные уравнения в порядок, чтобы получить оценки параметров модели методом наименьших квадратов. В этом случае легче использовать дискретное преобразование Фурье, чтобы обнаружить периодичности, сохранить только подмножество коэффициентов Фурье и инвертировать преобразование, чтобы получить установленные временные ряды.
Выполните спектральный анализ данных, чтобы выявить, какие частоты значительно способствуют изменчивости данных. Поскольку общее среднее значение сигнала приблизительно 9 000 и пропорционально преобразованию Фурье на 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.