Аccounting for autocorrelation in a linear regression problem on an example of analysis of atmospheric column NO2 content

Full Text

ВВЕДЕНИЕ

Одно из важных свойств многих геофизических процессов – память. Так, наличие дальней памяти значительно влияет на способность предсказания и/или детектирования долговременных изменений [1]. В регрессионном анализе временных рядов память может проявляться через автокорреляцию (сериальную корреляцию) остаточного ряда, и наличие автокорреляции влияет на погрешности регрессионных коэффициентов [2].

Известны несколько эффективных в практическом плане методов учета сериальной корреляции в оценках ошибок регрессионных коэффициентов. Они основаны на предположении о том, что остаточный ряд есть авторегрессиионный (АР) процесс первого порядка. В АР модели k-го порядка текущее значение процесса выражается в виде линейной комбинации k предыдущих значений процесса [3]. В АР процессе 1-го порядка текущее значение зависит только от одного, непосредственно предшествующего ему значения.

Наличие автокорреляции ведет к сокращению числа независимых значений, то есть, к уменьшению числа степеней свободы. В [4] оно выражено через АР коэффициент (равный коэффициенту автокорреляции при единичном сдвиге) и получена поправка к дисперсии остаточного ряда. Этот способ был, например, использован при оценке тренда температуры поверхности океана в [5]. В [6] предложена несколько иная аппроксимация дисперсии остаточного ряда с учетом автокорреляции, также выраженная через коэффициент АР процесса 1-го порядка. Данный метод был применен, например, в [7] при оценках трендов NO₂. В [8] предложена итерационная процедура учета сериальной корреляции в случае, когда остаточный ряд является АР процессом первого порядка (см. также [9]). Предполагается, что итерационный процесс сходится к глобальному минимуму суммы квадратов остаточного ряда. Эта процедура использовалась также в предположении, что остаточный ряд является АР процессом 2-го порядка [10]. Такой подход использовался и для учета сериальной корреляции в недавних оценках озонных трендов [11].

Автокорреляционная (АК) функция АР процесса 1-го порядка экспоненциально затухает [3]. Однако АК функция действительного остаточного ряда может значительно отличаться от экспоненты. Поэтому в рамках АР представления может потребоваться больший порядок авторегрессии для описания более сложного поведения остаточного ряда. Это означает, что память остаточного ряда проявляется на большем масштабе, чем предполагается АР процессом 1-го порядка, и это характерно для разных геофизических временных рядов.

Цель настоящей работы состоит в разработке метода учета сериальной корреляции на больших масштабах изменчивости остаточного ряда, который может быть представлен АР процессом высокого порядка. Эффективность метода демонстрируется на примере анализа данных 26-летних наземных спектрометрических измерений общего содержания NO₂ на Звенигородской научной станции ИФА им. А.М. Обухова РАН.

1. МЕТОД

Линейную регрессионную модель удобно записать в матричной форме [2]:

Y = Xβ + ε, (1)

где Y – вектор n наблюдений (анализируемый ряд) размера (n × 1), X – матрица предикторов (независимых переменных) размера (n × p), p – число предикторов, β – вектор неизвестных коэффициентов (искомых коэффициентов регрессии) размера (p × 1) , ε – вектор ошибок (остатков) длины n (остаточный ряд). Переопределенная система уравнений (1) (n > p) решается стандартным методом наименьших квадратов. Если остаточный ряд ε подчиняется гауссову распределению с нулевым средним и дисперсией σ², то стандартная ошибка искомых коэффициентов находится одновременно с решением системы [2].

1.1. Учет автокорреляции (сериальной корреляции)

Приведем, следуя [2], основные соотношения для решения задачи (1), то есть выражения для среднеквадратичных оценок коэффициентов регрессии и их стандартных ошибок. Пусть V(ε) есть ковариационная матрица остаточного ряда размера (n × n), такая что V(ε) = Vσ², где σ² – дисперсия некоррелированной компоненты ε. Если ошибки (члены остаточного ряда) ε некоррелированы, то V(ε) = Iσ², где I – единичная матрица, и матрица V(ε) есть просто диагональная дисперсионная матрица. Однако для геофизических данных типично, когда ошибки сериально коррелированны и, по крайней мере, некоторые недиагональные элементы матрицы V(ε) отличны от нуля.

Согласно [2], существует невырожденная симметричная матрица P, такая что

P^T P = P P = P² = V, (2)

где P^T – транспонированная матрица для матрицы P. Если остаточный ряд ε можно представить в виде АР процесса, то матрица V будет его АК матрицей [2]. Например, в случае АР процесса первого порядка элемент V_ij матрицы V есть ρ^|i–j|, где ρ – коэффициент автокорреляции при единичном запаздывании [3] (см. также приложение А в [9]).

Обозначим через e вектор некоррелированных ошибок, ассоциированных с вектором коррелированных ошибок ε. Тогда, согласно [2],

e = P^–1ε (3)

и искомый вектор b, элементы которого являются обобщенными среднеквадратичными оценками элементов вектора регрессионных коэффициентов β, есть:

b = (X^T V^–1 X)^–1 X^T V^–1 Y. (4)

При этом вектор квадрата стандартной ошибки δ_b² оценoк b состоит из диагональных элементов матрицы (X^T V^–1 X)^–1σ²:

δ_b² = diag(X^T V^–1 X)^–1σ². (5)

Формулы (4) и (5) описывают решение регрессионной задачи (1), обобщенное на случай сериальной корреляции. В отсутствие сериальной корреляции V^–1 = I и формулы (4), (5) принимают стандартный вид [2].

1.2. Метод решения

Как будет показано в разделе 4, низкий порядок АР представления остаточного ряда может быть недостаточным для учета долговременной корреляции. Вместо оценочных формул [4–7] и итерационной процедуры [8–11] мы предлагаем метод оценки регрессионных коэффициентов и их стандартных ошибок с использованием обобщенных формул (4), (5). Ключевым моментом является конструирование АК матрицы V АР процесса, аппроксимирующего поведение остаточного ряда. Известно, что АК функция АР процесса и коэффициенты авторегрессии связаны рекуррентными соотношениями, образующими систему уравнений Юла–Уокера (Yule–Walker) [3].

Рассмотрим два способа нахождения матрицы V. Пусть АР процесс имеет порядок k. АК функцию остаточного ряда для задержек от 1 до k рассчитываем обычным путем и, решая систему уравнений Юла–Уокера, по значениям АК коэффициентов находим k коэффициентов АР процесса. Затем по рекуррентным соотношениям рассчитываем остаток (хвост) АК функции (для задержек от k + 1 до n). В результате квадратная матрица V размера (n × n) полностью определена.

В более сложном, но более оптимальном способе используем алгоритм Бёрга (Burg) для расчета АР коэффициентов, который основан на принципе максимальной энтропии [12]. По значениям АР коэффициентов, путем решения системы уравнений Юла-Уокера относительно АК коэффициентов, находим k значений коэффициентов автокорреляции, соответствующих АР процессу, а затем рассчитываем прежним способом хвост АК функции.

Таким образом, регрессионная задача (1) решается в два шага. На первом шаге решение находится стандартным методом наименьших квадратов (в предположении V = I). Остаточный ряд представляется в виде АР процесса порядка k, который может значительно отличаться от 1. В алгоритме Бёрга k может достигать трети длины ряда. АР коэффициенты и значения соответствующей АК функции для задержек в интервале от 1 до k рассчитываются одним из описанных способов. Затем рассчитывается хвост АК функции и определяется матрица V.

На втором шаге в соответствии с (4) и (5) рассчитываются оценки регрессионных коэффициентов и их стандартные ошибки. При необходимости с помощью (3) может быть найден вектор некоррелированных ошибок e. Матрица P есть решение квадратного матричного уравнения (2) и находится с помощью преобразования матрицы V к жордановой форме.

В настоящей работе для анализа данных использован метод с применением алгоритма Бёрга.

2. ИСПОЛЬЗОВАННЫЕ ДАННЫЕ И РЕГРЕССИОННАЯ МОДЕЛЬ

Для демонстрации предложенного метода воспользуемся данными измерений общего содержания NO₂ на Звенигородской научной станции в 1990–2016 гг. Станция расположена в сельской местности в 50 км западнее Москвы и входит в состав Сети по обнаружению изменений состава атмосферы (Network for the Detection of Atmospheric Composition Change – NDACC). Измерения выполняются в утренние и вечерние сумерки по рассеянной из зенита солнечной радиации в видимой области спектра. Метод наблюдений и восстановления содержания NO₂ описан в [13–15], а данные находятся в свободном доступе по адресу ftp://ftp.cpc.ncep.noaa.gov/ndacc/station/zvenigor/. Анализируемое общее содержание NO₂ есть интегральное содержание NO₂ в вертикальном столбе атмосферы выше ее пограничного слоя. В настоящей работе использованы данные вечерних измерений. Результаты измерений подверглись статистической отбраковке по способу [16], и по ним рассчитаны среднемесячные значения содержания NO₂. Они приведены на рис. 1а.

В модель множественной линейной регрессии для NO₂ включены следующие предикторы: линейный член, зональная скорость экваториального ветра на изобарической поверхности 20 гПа и среднее зональных скоростей на поверхностях 40 и 50 гПа (http://www.geo.fu-berlin.de/en/met/ag/strat/produkte/qbo/index.html) – в качестве индексов квазидвухлетней цикличности (КДЦ), индекс Ниньо 3.4 (https://www.esrl.noaa.gov/psd/data/correlation/nina34.data) и смещенный во времени индекс Ниньо 3.4 – для учета эффектов Эль-Ниньо – Южного колебания (ЭНЮК), индекс Североатлантического колебания (САК) (https://crudata.uea.ac.uk/cru/data/nao/nao.dat), индекс солнечной активности F10.7 (https://www.ngdc.noaa.gov/stp/space-weather/solar-data/solar-features/solar-radio/), среднее по северному полушарию значение оптической толщи стратосферного аэрозоля на длине волны 550 нм (https://data.giss.nasa.gov/modelforce/strataer/) – для учета воздействия продуктов извержения вулкана Пинатубо на NO₂. Эти предикторы представляют в регрессионной модели основные важные процессы, оказывающие влияние на изменчивость стратосферных примесей. В различных комбинациях они использовались при анализе трендов NO₂ [7, 16–21].

Два индекса КДЦ использованы для учета задержки отклика NO₂ на воздействие экваториальной КДЦ [22–23]. Вариации скоростей ветра на изобарической поверхности 20 гПа и в слое 40–50 гПа в экваториальной стратосфере в высокой степени ортогональны между собой (коэффициент корреляции КК ~ –0.01). Использование двух ортогональных индексов КДЦ также типично при анализе озонных трендов [10–11].

Аналогичная идея использована и для учета эффекта ЭНЮК. Значения второго индекса запаздывают относительно значений оригинального индекса Ниньо 3.4 на 14 месяцев, так что два индекса ортогональны между собой (КК ~ 0.01).

Стоит отметить, что не все предикторы являются попарно независимыми. В наибольшей степени это относится к аэрозольной оптической толще (АОТ). Форму ее временной зависимости с обратным знаком можно представить по штриховой кривой на рис. 1б. Максимум толщи, обусловленный продуктами извержения вулкана Пинатубо, приходится на начало 1992 г. Поэтому изменения АОТ в период 1990–2016 гг. имеют трендовую составляющую. Кроме того, максимум АОТ пришелся на окончание фазы максимума солнечного цикла, и уменьшение АОТ и солнечного индекса F10.7 в последующие пять лет происходили синхронно, за счет чего между этими предикторами есть корреляционная связь. Обусловленный этим КК ~ 0.1 для периода 1990–2016 гг. Корреляционная связь с АОТ может возникнуть при наложении максимума АОТ на экстремум другого предиктора, как это случилось с индексом Ниньо 3.4 (КК = –0. 27). Однако исключать АОТ из регрессионной модели не целесообразно, так как извержение вулкана Пинатубо привело к значительному влиянию на стратосферное содержание NO₂ [7, 13, 16–18, 20, 24].

Наличие корреляции между предикторами приводит к корреляции между оценками регрессионных коэффициентов (4), и значения парных КК определяются недиагональными элементами матрицы (X^T V^–1 X)^–1 [2]. Корреляция между предикторами отражается и на ошибках регрессионных коэффициентов (см. (5)).

Все искомые коэффициенты в регрессионной модели представлены в виде разложения в ряд Фурье по парам синус-косинус, соответствующим годовой и полугодовой гармоникам, с целью учета годового хода NO₂, сезонной зависимости трендов NO₂ и сезонной зависимости воздействия других предикторов на NO₂ (см. также [10–11]). Годовой ход NO₂ на Звенигородской станции в основном обусловлен годовой и полугодовой гармониками [13].

3. РЕЗУЛЬТАТЫ АПРОБАЦИИ МЕТОДА НА ПРИМЕРЕ NO₂

Рисунок 1а показывает, что содержание NO₂, предсказываемое регрессионной моделью, хорошо соответствует наблюдаемому содержанию NO₂ (КК = 0.99). Это связано в первую очередь с большой амплитудой годового хода NO₂, которая намного превышает амплитуды межмесячных и межгодовых вариаций. С исключенным годовым ходом из наблюдаемых и модельных данных КК = 0.63.

Как показывает штриховая кривая на рис. 1б, значительное уменьшение содержания NO₂ в 1991–1992 гг. есть результат влияния на NO₂ продуктов извержения вулкана Пинатубо (в июне 1991 г.). Оно обусловлено гетерогенными химическими процессами на поверхности стратосферного сульфатного аэрозоля [25]. Однако стоит отметить, что индекс АОТ недооценивает вулканический эффект в 1993 г. (ср. две кривые на рис. 1а и сплошную черную кривую на рис. 1б). Согласно наблюдениям, восстановление NO₂ происходило медленнее, чем предсказывается моделью с АОТ в качестве предиктора.

Коррелированная ошибка (остаток) ε и некоррелированная ошибка e, соответствующая порядку авторегрессии 50, приведены на рис. 1б, а на рис. 2 сравниваются спектры ошибок (рассчитаны с помощью метода максимальной энтропии [12]). Спектр коррелированной ошибки указывает на наличие в остаточном ряде ε вариаций с периодом около 80 мес. Спектр имеет общий наклон; черная штриховая прямая – степенная аппроксимация с показателем степени 0.7. В отличие от него спектр некоррелированной ошибки следует степенному закону с показателем, близким к нулю, что является теоретической величиной для спектра белого шума.

На рис. 3 показаны АК функции АР процессов, аппроксимирующих временной ряд коррелированной ошибки при различных значениях порядка авторегрессии k. АК функция, соответствующая значению k = 90 (сплошная серая кривая), демонстрирует колебания с периодом около 80 мес, что согласуется с периодом, проявившимся в спектре на рис. 2. АК функция для k = 1 экспоненциально затухает, а в случае k = 10 затухание проявляется на задержках более 10 мес. При значении k = 30 АК функция содержит колебания, но их период пока меньше 80 мес. Значение k = 50 достаточно, чтобы разрешить 80-месячные вариации, и соответствующая АК функция, совпадая в диапазоне k от 1 до 50 мес с функцией для k = 90, в целом следует ей и при бóльших задержках.

Для демонстрации влияния порядка авторегрессии на оценки коэффициентов регрессии, на рис. 4 приведены годовые оценки линейного тренда NO₂ и их 95%-е доверительные интервалы в зависимости от порядка k АР процесса. Значение k = 0 соответствует случаю без учета сериальной корреляции. Рисунок 4 показывает, что величина k влияет и на значения, и на ошибки регрессионных коэффициентов. Оценки тренда при разных значениях k изменяются примерно на треть его величины; однако различия находятся в пределах 95%-х доверительных интервалов. Ошибка оценки тренда значительно, более чем вдвое, увеличивается с увеличением k при небольших значениях k, достигает максимума при k = 12÷15 и постепенно уменьшается, пока k не достигает 50. При k > 50 оценки тренда и их ошибки остаются примерно постоянными.

На рис. 4 серым цветом показана также оценка тренда и доверительный интервал, рассчитанные по методу [8] (итерационная процедура). Крестиком показан тренд и доверительный интервал, оцененный по методу [6]. Эти оценки близки между собой, а значения трендов отличаются от оценки тренда, полученного по нашему методу при k = 1. Эти значения трендов практически (для метода [6] в точности) совпадают со значением тренда при k = 0, то есть в отсутствие сериальной корреляции.

На рис. 5 приведены оценки линейных трендов NO₂ и величин изменений содержания NO₂, обусловленных различными предикторами, полученные при k = 50. На рисунках нанесены годовые, месячные и сезонные значения. Процентные значения рассчитаны относительно многолетнего среднегодового значения содержания NO₂. Линейный тренд на рис. 5а отражает изменение содержания NO₂ за 10 лет. Эффект 11-летнего солнечного цикла (рис. 5б) соответствует изменению содержания NO₂ при изменении уровня солнечной активности в 11-летнем солнечном цикле от минимума к максимуму. Вулканический эффект на рис. 5в означает изменение содержания NO₂, соответствующее максимуму АОТ. Эффект КДЦ на рис. 5г есть изменение содержания NO₂ при изменении экваториальной зональной скорости ветра от максимальной восточной до максимальной западной фазы. Эффекты САК и ЭНЮК (рис. 5д, е) – это изменения содержания NO₂, обусловленные изменениями индексов САК и Ниньо 3.4 соответственно, на величины их среднеквадратичных отклонений. Эффект ЭНЮК приведен в виде отклика NO₂ на изменение только одного индекса Ниньо 3.4, так как отклик на другой, ортогональный ему индекс статистически незначим на уровне 95%.

Согласно рис. 5а, линейный тренд NO₂ сезонно зависим с максимальными по модулю значениями весной (до –6% за 10 лет) и осенью. Годовой тренд составляет около –3% за 10 лет и также статистически значим.

В то время как годовая оценка эффекта солнечной активности в NO₂ невелика и статистически незначима, значительный эффект получен для периода с конца зимы по середину весны (рис. 5б). Так, содержание NO₂ в марте во время фазы минимума солнечной активности на 8% больше, чем во время фазы максимума. Относительно среднемесячных значений содержания NO₂ эффект солнечного цикла в феврале и марте составляет около –10%.

Максимальный абсолютный эффект извержения вулкана Пинатубо, около –30%, получен для летнего сезона (рис. 5в). Однако если изменения содержания NO₂ нормировать на среднемесячные его значения, то вулканический эффект в меньшей степени зависит от сезона и варьирует вблизи значения –20%.

Эффект КДЦ в NO₂ выражен преимущественно в весеннее-летний период (рис. 5г). Его величина составляет летом около 8% относительно среднегодового содержания NO₂. Разный сезонный ход двух кривых на рис. 5г указывает на то, что фазовый сдвиг эффекта КДЦ в NO₂ относительно экваториальной КДЦ изменяется с сезоном. Например, эффект в NO₂ в июне совпадает по фазе со скоростью ветра на 20 гПа. Эффект КДЦ в апреле запаздывает относительно КДЦ экваториальной скорости ветра на уровне 20 гПа на одну восьмую периода КДЦ (~3.5 мес), поскольку апрельские значения двух кривых на рис. 5г почти одинаковы, а фаза КДЦ экваториального ветра распространяется сверху вниз. Летние значения содержания NO₂ во время западной (положительной) фазы экваториального ветра на уровне 20 гПа больше, чем во время восточной фазы. Весенние значения содержания NO₂ во время западной фазы КДЦ во всем слое 20–50 гПа в целом больше, чем во время восточной фазы.

Небольшой, отрицательный по величине, статистически значимый эффект САК в NO₂ получен для весеннего и, на пределе статистической значимости, осеннего периодов (рис. 5д). Изменения содержания NO₂ весной оцениваются в 2–3% при изменении индекса САК на величину его среднеквадратичного отклонения.

Отрицательный эффект ~3% получен для эффекта ЭНЮК в конце зимы – начале весны, а эффект противоположного знака, но на пределе статистической значимости, проявился в июне (рис. 5е).

Следует отметить, что оценки линейных трендов NO₂, эффекта 11-летнего солнечного цикла и вулканического эффекта в NO₂ на Звенигородской станции, полученные методом множественной линейной регрессии, были представлены и в более ранних работах [16, 18, 26]. Знаки оценок из настоящей работы и из указанных работ совпадают. Годовая оценка тренда NO₂, полученная с помощью нового метода, близка к наиболее поздней оценке в [26], хотя и несколько меньше ее по модулю, но в два раза меньше (по модулю) более ранних оценок в [16, 18]. Вулканический эффект по результатам настоящей работы близок к эффектам, полученным в [16, 18, 24] по данным вечерних измерений, но в полтора раза меньше приведенных в [27, 28] ранних оценок, полученных без использования множественной регрессии. Для эффекта солнечного цикла в вечерних данных об NO₂ в [24] получено статистически значимое отрицательное значение, почти двукратно превышающее по модулю статистически незначимую годовую оценку, полученную в настоящей работе. Причины отмеченных количественных расхождений могут быть разными. Имеют значение как длина анализируемого ряда, так и различия использованных регрессионных моделей, в частности, выбор предиктора для описания влияния вулканического аэрозоля.

4. ЗАКЛЮЧЕНИЕ

В работе предложен метод учета сериальной корреляции данных в задаче линейной регрессии. Одно из достоинств метода – то, что он позволяет учесть автокорреляцию остаточного ряда, которая может проявляться на больших масштабах. Эффективность метода продемонстрирована на основе анализа данных 26-летних измерений общего содержания NO₂ на Звенигородской научной станции с применением модели множественной линейной регрессии.

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

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

Сначала ошибка увеличивается с увеличением k. Затем она может достичь максимума и после этого начать уменьшаться. В случае с NO₂ на Звенигородской станции ошибка более чем удваивается в своем максимуме по сравнению с ее начальным значением.

Уменьшение ошибки может прекратиться при больших значениях k, если k достигнет значения, при котором АР процесс этого порядка обеспечит воспроизведение важных особенностей автокорреляционной функции остаточного ряда. Однако если ошибка не достигает максимума в диапазоне допустимых значений k, то следует полагать, что использованная регрессионная модель не вполне адекватно описывает поведение наблюдаемой переменной и должна быть усовершенствована (см. также [2]).

В случае с использованными данными об NO₂ максимум ошибки пришелся на k~10, а затем она уменьшалась по достижении k значения ~50. В диапазоне k от 50 до 90 регрессионные коэффициенты и их ошибки сохраняли примерно постоянные значения.

Модель множественной линейной регрессии для данных наблюдений общего содержания NO₂ включает сезонно зависимый годовой тренд и сезонно зависимые эффекты воздействия на NO₂ 11-летнего солнечного цикла, экваториальной КДЦ, САК, ЭНЮК и продуктов извержения вулкана Пинатубо. Годовые и сезонно зависимые оценки вклада этих эффектов в долговременную и межгодовую изменчивость NO₂ получены с учетом автокорреляции ошибки (остаточного ряда) до временного масштаба 50 мес.

Представленный в работе метод применим к анализу самых разных данных. Количественные характеристики зависимости регрессионных коэффициентов и их ошибок от порядка АР процесса для них будут своими.

Источник финансирования. Работа выполнена при финансовой поддержке РФФИ (проект № 16-05-00663) и программ Российской академии наук. Данные измерений содержания NO₂ архивированы NDACC. Данные об индексе солнечной активности F10.7 представлены для свободного доступа Национальным геофизическим центром NOAA (NOAA National Geophysical Data Center). Данные о зональной скорости экваториального стратосферного ветра подготовлены Свободным берлинским университетом (Freie Universität Berlin). Индекс Ниньо 3.4 представлен Лабораторией исследований земной системы NOAA (NOAA-ESRL Physical Sciences Division). Данные об индексе САК представлены Климатическим отделом Университета Восточной Англии (Climatic Research Unit, University ofEast Anglia). Данные об аэрозольной оптической толще подготовлены Годдардовским институтом космических исследований НАСА (NASA Goddard Institute for Space Studies).

About the authors

A. N. Gruzdev

Author for correspondence.
Email: a.n.gruzdev@mail.ru
Russian Federation, Pyzhevsky per., 3, Moscow, 119017

References

Supplementary files