Линейная регрессия частотного диапазона

В этом примере показано, как использовать дискретное преобразование Фурье для создания линейной регрессионой модели для временных рядов. Временные ряды, используемые в этом примере, являются ежемесячным количеством случайных смертей в Соединенных Штатах с 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 месяцев). Периодический характер данных предполагает, что соответствующая модель может быть

где $$\mu$$ - общее среднее, $$N$$ - длина временных рядов, и $$\epsilon (n)$$ является белой шумовой последовательностью независимых и одинаково распределенных (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.

См. также

fft | ifft | xcorr