4.1: Вступ до спектрального аналізу

Last updated
Save as PDF

Page ID: 97222

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)

Багато часових рядів, обговорюваних у попередніх розділах, відображали сильні періодичні компоненти: номери сонячних плям Приклад 1.1.1, кількість захоплених рисі Прикладу 1.1.2 та австралійські дані про продаж вина Приклад 1.4.1. Часто існує очевидний вибір на період\(d\) такої циклічної частини, як річна картина продажів вина. Враховуючи\(d\), можна було б продовжити, видаляючи сезонні ефекти, як у розділі 1.4. У перших двох прикладах, однак, дещо складніше визначити точне значення\(d\). У цьому розділі обговорюється загальний метод для розгляду періодичних компонентів часового ряду. Щоб ускладнити справу, зазвичай трапляється, що в часовому ряду одночасно присутні кілька циклічних закономірностей. Як приклад нагадаємо дані південного індексу коливань (SOI), які демонструють як річну картину, так і так звану модель El Ni\(\tilde{n}\) o.

Функції синуса і косинуса є прототипами періодичних функцій. Вони будуть використані тут для опису циклічної поведінки в часових рядах. Перш ніж це зробити, цикл визначається як один повний період синуса або косинуса функції протягом часового інтервалу довжини\(2\pi\). Визначте також частоту

\[ \omega=\dfrac 1d \nonumber \]

як кількість циклів за спостереження, де\(d\) позначається період часового ряду (тобто кількість спостережень в циклі). Для щомісячних спостережень з річним періодом,\(d=12\) а значить і\(\omega=1/12=0.083\) циклами за спостереження. Тепер перегляньте процес

\[ X_t=R\sin(2\pi\omega t+\varphi) \nonumber \]

як введено в прикладі 1.2.2, використовуючи конвенцію\(\lambda=2\pi\omega\). Щоб включити випадковість в цей процес, виберіть амплітуду\(R\) і фазу,\(\varphi\) щоб бути випадковими величинами. Еквівалентне представлення цього процесу дається

\[ X_t=A\cos(2\pi\omega t)+B\sin(2\pi\omega t), \nonumber \]

з\(A=R\sin(\varphi)\) і,\(B=R\cos(\varphi)\) як правило, незалежні стандартні нормальні варіації. Потім,\(R^2=A^2+B^2\) є\(\chi\) -квадрат випадкової величини з 2 ступенями свободи і\(\varphi=\tan^{-1}(B/A)\) рівномірно розподіляється по\((-\pi,\pi]\). Причому,\(R\) і\(\varphi\) є незалежними. Вибравши зараз значення\(\omega\) однієї конкретної періодичності можна описати. Для розміщення більше одного, здається природним розглянути суміші цих періодичних серій з декількома частотами і амплітудами:

\[ X_t=\sum_{j=1}^m \big[A_j\cos(2\pi\omega_jt)+B_j\sin(2\pi\omega_jt)\big], \qquad t\in\mathbb{Z}, \nonumber \]

де\(A_1,\ldots,A_m\) і\(B_1,\ldots,B_m\) є незалежними випадковими величинами з нульовим середнім і\(\sigma_1^2,\ldots,\sigma_m^2\) дисперсіями, і\(\omega_1,\ldots,\omega_m\) є різними частотами. Можна показати, що\((X_t\colon t\in\mathbb{Z})\) є слабо стаціонарним процесом з лаг- ч ACVF

\[ \gamma(h)=\sum_{j=1}^m\sigma_j^2\cos(2\pi\omega_j h),\qquad h\in\mathbb{Z}. \nonumber \]

Останній результат дає, зокрема, що\(\gamma(0)=\sigma_1^2+\ldots+\sigma_m^2\). Дисперсія\(X_t\) - це, отже, сума дисперсій складових.

Приклад 4.1.1. Нехай\(m=2\) і вибирайте\(A_1=B_1=1\),\(A_2=B_2=4\) щоб бути постійним так само, як\(\omega_1=1/12\) і і\(\omega_2=1/6\). Це означає, що

\[ X_t=X_t^{(1)}+X_t^{(2)} =\big[\cos(2\pi t/12)+\sin(2\pi t/12)\big]+\big[4\cos(2\pi t/6)+4\sin(2\pi t/6)\big] \nonumber \]

це сума двох періодичних складових, з яких одна демонструє річний цикл, а інша - цикл на півроку. Для всіх задіяних процесів реалізації\(n=48\) спостережень (дані за 4 роки) відображені на малюнку 4.1. Також показаний графік четвертого часового ряду, який містить\(X_t\) спотворений стандартним нормальним незалежним шумом,\(\tilde X_t\). Відповідний код R такий:

> t = 1:
48> x1 = cos (2* пі* т/12) +гріх (2* пі* т/12)
> x2 = 4* cos (2* пі* т/6) +4* гріх (2* пі* т/6)
> х = х1+х2
> тильдекс = х+норма (48)

Зверніть увагу, що квадрат амплітуди\(X_t^{(1)}\) дорівнює\(1^2+1^2=2\). Таким чином, максимальне і\(X_t^{(1)}\) мінімальне значення є\(\pm\sqrt{2}\). Аналогічно отримуємо і\(\pm\sqrt{32}\) для другого компонента.

Для статистика зараз важливо розробити інструменти для відновлення періодичності даних. Гілка статистики, що займається цією проблемою, називається спектральним аналізом. Стандартний метод в цій області заснований на періодограмі, яка введена зараз. Припустимо на той момент, що параметр частоти\(\omega_1=1/12\) в прикладі 4.1.1 відомий. Для отримання оцінок\(A_1\) і\(B_1\), можна спробувати запустити регресію за допомогою пояснювальних змінних\(Y_{t,1}=\cos(2\pi t/12)\) або\(Y_{t,2}=\sin(2\pi t/12)\) обчислити оцінювачі найменших квадратів

\ почати {вирівнювати*}
\ капелюх A_1=&\ dfrac {\ сума_ {t=1} ^nx_ty_ {t,1}} {\ сума {t=1} ^nY_ {t,1} ^2} =\ dfrac 2n\ sum_ {t=1} ^nx_t\ cos (2\ пі т/12),\\ [.2см]
\ капелюх B_1=&\ dfrac {\ сума_ {t=1} ^nx_ty_ {t,2}} {\ sum_ {t=1} ^nY_ {t,2} ^2} =\ dfrac 2n\ sum_ {t=1} ^nx_t\ sin (2\ пі т/12).
\ end {вирівнювати*}

Малюнок 4.1: Графіки часових _{рядів (X t} ⁽¹⁾), (X _t ⁽²⁾), (X _t) та\(\tilde X_t\)

Оскільки, в цілому, задіяні частоти не будуть відомі статистику до аналізу даних, вищесказане пропонує вибрати ряд потенційних\(\omega's, say j/n for \(j=1,\ldots,n/2\) і запустити тривалу регресію форми

\[X_t=\sum_{j=0}^{n/2}\big[A_j\cos(2\pi jt/n)+B_j\sin(2\pi jt/n)\big]. \tag{4.1.1} \]

Це призводить до найменших квадратів оцінок\(\hat A_j\) і\(\hat B_j\) з яких слід вибирати «значущі». Зверніть увагу, що регресія в 4.1.1 є ідеальним, тому що існує стільки невідомих, скільки змінних! Зауважте також, що

\[ P(j/n)=\hat A_j^2+\hat B_j^2 \nonumber \]

є по суті (аж до нормалізації) оцінювачем кореляції між часовим рядом\(X_t\) і відповідною сумою періодичних косинусних і синусоїдних функцій на частоті\(j/n\). Колекція всіх\(P(j/n)\)\(j=1,\ldots,n/2\), називається масштабована періодограма. Його можна швидко обчислити за допомогою алгоритму, відомого як швидке перетворення Фур'є (БПФ), який, в свою чергу, заснований на дискретному перетворенні Фур'є (DFT)

\[ d(j/n)=\dfrac{1}{\sqrt{n}}\sum_{t=1}^nX_t\exp(-2\pi ijt/n). \nonumber \]

Частоти\(j/n\) називаються Фур'є або основними частотами. Так як\(\exp(-ix)=\cos(x)-i\sin(x)\) і\(|z|^2=z\bar{z}=(a+ib)(a-ib)=a^2+b^2\) для будь-якого комплексного числа випливає\(z=a+ib\), що

\[ I(j/n)=|d(j/n)|^2=\dfrac 1n\left(\sum_{t=1}^nX_t\cos(2\pi jt/n)\right)^2+\dfrac 1n\left(\sum_{t=1}^nX_t\sin(2\pi jt/n)\right)^2. \nonumber \]

Кількість\(I(j/n)\) позначається як періодограма. Відразу випливає, що періодограма і масштабована періодограма пов'язані через ідентичність\(4I(j/n)=nP(j/n)\).

Приклад 4.1.2. Використовуючи вирази та позначення Прикладу 4.1.1, періодограма та масштабована періодограма обчислюються в R наступним чином:

> t = 1:48
> l = abs (fft (x) /sqrt (48)) ^
2> P = 4*
I/48> f = 0:24
/48>ділянка (f, P [1:25] , тип = «л»)
> аблайн (v=1/12)
> аблайн (v=1/6)

Відповідну (масштабовану) періодограмму для\((\tilde X_t)\) можна отримати аналогічним чином. Масштабовані періодограми показані на лівій і середній панелі малюнка 4.2. На правій панелі відображається масштабована періодограма іншої версії,\((\tilde X_t)\) в якій стандартний нормальний шум був замінений звичайним шумом з дисперсією 9. З цих графіків видно, що періодичність шести місяців добре видно на графіках (див. Пунктирні вертикальні лінії при x=1/6. Менш виражений річний цикл (вертикальна лінія при x = 1/12 все ще видно на перших двох масштабованих періодограмах, але втрачається, якщо дисперсія шуму збільшується, як на правій графіці. Однак зауважте, що шкала y відрізняється для всіх трьох графіків.

Малюнок 4.2: Масштабовані періодограми (\(X_t\)),\((\tilde X_t^{(1)})\),\((\tilde X_t^{(2)})\)

В ідеальній ситуації, коли ми спостерігаємо періодичну складову без додаткового забруднення шумом, ми можемо, крім того, зрозуміти, чому періодограма може бути корисною для виявлення дисперсійного розкладання зверху. Ми показали в рядках, що передують Приклад 4.1.1, що квадратні амплітуди\(X_{t}^{(1)}\) і\(X_t^{(2)}\) складають 2 і 32 відповідно. Ці значення легко зчитуються з масштабованої періодограми на лівій панелі малюнка 4.2. Забруднення шумом змінює ці значення.

У наступному розділі встановлено, що підхід в часовій області (заснований на властивостях ACVF, тобто регресії за минулими значеннями часового ряду) та підхід частотної області (з використанням періодичного функціонального підходу через фундаментальні частоти, тобто регресія на синусоїдних і косинусних функціях) еквівалент. Деякі подробиці наведено про спектральну щільність (аналог популяції періодограми) і про властивості самої періодограми.