CC BY
154
УДК 597.553.2-154.343 : 57.087.1
ПРОГНОЗИРОВАНИЕ ПОДХОДОВ ЛОСОСЕВЫХ (НА ПРИМЕРЕ КИЖУЧА ЗАПАДНОЙ КАМЧАТКИ) С ИСПОЛЬЗОВАНИЕМ МОДЕЛЕЙ ЭКСТРАПОЛЯЦИИ ВРЕМЕННЫХ РЯДОВ И МОДЕЛЕЙ «ЗАПАС-ПОПОЛНЕНИЕ»
М.Г. Фельдман, Е.А. Шевляков, Ж.Х. Зорбиди
Ст. н. с., зав. лаб., вед. н. с., Камчатский научно-исследовательский институт рыбного хозяйства и океанографии
683000 Петропавловск-Камчатский, Набережная, 18 Тел., факс: (4152) 41-27-01, (4152) 42-07-74 E-mail: feldman.m.g.@kamniro.ru
ТИХООКЕАНСКИЕ ЛОСОСИ, КИЖУЧ, ПРОГНОЗИРОВАНИЕ, ВРЕМЕННЫЕ РЯДЫ, АВТОРЕГРЕССИЯ, «ЗАПАС-ПОПОЛНЕНИЕ», ARIMA, МОДЕЛЬ РАСПРЕДЕЛЕННОГО ЛАГА
В данной работе рассматривается проблема прогнозирования подходов тихоокеанских лососей. В качестве модельной популяции взято стадо кижуча Oncorhynchus kisutch Западной Камчатки. В регрессионном анализе использовались и сравнивались модели двух классов: модели временных рядов и модели зависимости пополнения популяции от нерестового запаса. Предлагается представлять прогноз либо в виде произведения плотностей вероятности прогнозов по различным моделям, либо как средневзвешенный результат плотностных прогнозов из нескольких регрессионных моделей.
FORECASTING THE RUNS OF PACIFIC SALMON (FOR EXAMPLE - A COHO SALMON OF THE WESTERN KAMCHATKA) ON THE BASE OF “STOCK-RECRUITMENT” MODELS AND TIME SERIES EXTRAPOLATION SIMULATIONS M.G. Feldman, E.A. Shevlyakov, Zh.Kh. Zorbidi
Senior scientist, head of dep., leading scientist, Kamchatka Research Institute of Fisheries and Oceanography 683000 Petropavlovsk-Kamchatsky,Naberedzhnaya, 18 Tel., fax: (4152) 41-27-01, (4152) 42-07-74 E-mail: feldman.m.g.@kamniro.ru
PACIFIC SALMON, COHO SALMON, FORECASTING, TIME SERIES, AUTOREGRESSION, “STOCK-RECRUITMENT”, ARIMA, DISTRIBUTED LAG
The paper provides analysis of the issues of forecasting runs of Pacific salmon. The stock used as a basis in simulation was the West Kamchatkan stock of coho salmon Oncorhynchus kisutch. Simulations of two classes were compared and used in regression analysis: of time series and of the relationship between the spawning stock and the recruitment. It is suggested to perform the forecast either as a product of probability density of forecasts of different models or as an average weighted result of density-based forecasts of several regression models.
Прогнозирование динамики численности рыб, в промысловой F и естественной смертности M, частности тихоокеанских лососей, представляет было получено: собой важную проблему рыбного хозяйства и име- Nt+1 = Nte~Zt.
ет достаточно долгую историю. Основоположни- Известно, что один из переводов данной рабо-
ком исследований популяционной динамики в ты получили Р. Бивертон и С. Холт (по данным
рыбном хозяйстве по праву считается Ф.И. Бара- Anderson, 2002) в 1943 г., а другой был сделан
нов. Его также называют «дедушкой рыбных по- У Рикером в 1945 г. (по данным Quinn, 2003).
пуляционных исследований» (Quinn, 2003). Наи- В дальнейшем У. Рикер (Ricker, 1954) предло-
более известной и важной его работой является жил решать задачу определения численности бу-
статья «К вопросу о биологических основаниях дущего подхода путем расчета пополнения исходя
рыбного хозяйства» (1918), в которой автор впер- из родительского нерестового запаса и представил
вые применяет дифференциальное счисление для первую математическую модель зависимости по-
моделирования динамики численности рыбного полнения от родительского запаса. В данном слу-
запаса. Решая уравнение: чае исходное дифференциальное уравнение для
dN / dt = -ZN , численности поколения имело вид:
где N — численность рыб, а Z=F+M — коэффици- dN / N = —(q + pS)dt, (1)
ент смертности, складываемый из коэффициентов где S — исходный нерестовый запас, q иp — неза-
висимая и зависимая от S компоненты смертности. Учитывая, что исходная численность поколения равна численности родительского запаса, умноженной на плодовитость — N=fS, решая уравнение (1), получаем:
R = aSe-bS, (2)
где S — нерестовый запас, R — ожидаемое пополнение, a и b — коэффициенты депенсационной и компенсационной смертности соответственно.
Графиком модели (2) является куполообразная кривая, показывающая снижение пополнения при высоких значениях нерестового запаса.
Немного позже, в 1957 г., Р. Бивертон и С. Холт (Бивертон, Холт, 1969) предложили альтернативную модель, где коэффициент смертности Z зависит не от родительского запаса S, как в (1), а от текущего значения численности N:
dN / N = -(q + pN)dt (3)
Решая уравнение (3) относительно численности пополнения для любого момента времени t, получим:
R = aS/(1 + (a / b)S) (4)
При этом a — максимум R/S при низком запасе, b — максимум R при высоком запасе (Хил-борн, Уолтерс, 2001). В отличие от кривой Рикера, кривая Бивертона-Холта представляет собой непрерывно возрастающую к горизонтальной асимптоте функцию.
В дальнейшем происходит развитие теории: Д. Кушинг (Cushing, 1973), Дж. Полик (Paulik, 1973) и Р. Деризо (Deriso, 1980) предлагают свои модели и вводят в них степенные параметры.
В 1982 г. Дж. Шепард (Sheperd, 1982), обобщая ранее полученные результаты, приходит к выводу
о единой структуре всех предложенных моделей: R = aSf (S / K'), где f(S/K') — функция убывания, зависящая от запаса, и К' — пороговая биомасса, при превышении которой в запасе начинают преобладать зависящие от плотности процессы. В качестве такой функции Дж. Шепард, стремясь получить более гибкую кривую, использовал степенной параметр c и предложил:
f(S/K)=1/(1+S/K')c, соответственно пополнение будет равно:
R = aS/(1 + S / K ')c (5)
При c<1 уравнение (5) аналогично модели Кушинга, при c=1 принимает вид модели Бивертона-Холта, а при с>1 модель становится аналогичной зависимости Рикера.
Более подробно с историей развития моделирования согласно теории «запас-пополнение» можно ознакомиться в работах В.К. Бабаяна (1988) и Е.А. Криксунова (1995).
Несмотря на зачастую довольно приблизительные результаты, получаемые по моделям «запас-пополнение», теория, с одной стороны, позволила установить безопасные пределы промыслового изъятия, с другой стороны, дала мощный толчок рыбохозяйственной статистике — сбору первичного материала, обязательному определению возрастного состава рыб в уловах, определению плодовитости, оценке нерестового запаса и т. п.
Таким образом, при решении задачи определения численности будущего подхода не обойтись без моделей, основанных на теории «запас-пополнение». В данной работе мы хотим сравнить результаты, полученные с помощью моделей этого класса, с другими подходами. В качестве таковых были выбраны модели, относящиеся к методам экстраполяции временного ряда и допустимые в прогнозировании динамики рыбных популяций (Бабаян, 2000), а именно модель авторегрессии и проинтегрированного скользящего среднего (autoregression and integrated moving average—ARIMA) и модель распределенного лага (distributed lag model — DLM).
Первая из них была предложена Дж. Боксом и Г. Дженкинсом (1974) и позволяет выделить в исследуемом временном ряду регулярные компоненты (тренд и циклические компоненты). Данный метод чрезвычайно популярен в эконометрических исследованиях, и практика подтвердила его мощность и гибкость (Hoff, 1983; Pankratz, 1983; Vandaele, 1983). Модель можно представить в следующем виде:
Д "X, = с + £(,, Д X, - +&,<;„ + е,, (6)
i=1 j =1
где Xt — временной нестационарный ряд, Ad — оператор разности временного ряда порядка d (последовательное взятие d раз разностей первого порядка — сначала от временного ряда, затем от полученных разностей первого порядка, затем от второго порядка и т. д.), с — константа, a. и b —
z J
параметры авторегрессии и скользящего среднего, st — ошибки модели.
Модель распределенного лага можно представить в общем случае как:
J, = С + b0+ b1 -1 +... + bnxn-1 + , (7)
где yt — зависимая исследуемая переменная, xt —
независимая объясняющая переменная, с — константа, £{ — ошибки модели.
Применительно к популяциям кижуча Западной Камчатки, проблема прогнозирования подходов данного вида рыб встает достаточно остро.
Сравнивая динамику заполнения нерестилищ производителями кижуча с его общим подходом (рис. 1), можно отметить парадоксальную ситуацию: начиная с середины 2000-х годов, при снижающемся пропуске производителей кижуча на нерестилища, его подходы начали расти, также возрастает и интенсивность промысла. Следовательно, увеличивается и эффективность воспроизводства этой рыбы, причем многократно, особенно если учесть при этом и значительное нелегальное изъятие (Запорожец и др., 2007). Действительно, согласно теоретическим моделям зависимости пополнения от запаса, эффективность воспроизводства (показатель отношения числа рекрутов на единицу производителя) растет при снижении плотности запаса. Однако при этом общее
количество потомков снижается. Таким образом, объяснить значительный рост эффективности воспроизводства можно с помощью одновременного влияния двух факторов: улучшения условий внешней среды, способствующих более высокой выживаемости в посткатадромный период жизни (Шевляков и др., 2013), и недоучета производителей на нерестилищах.
В пользу первого фактора косвенно говорит, например, одновременный рост подходов других видов тихоокеанских лососей, начиная с 2006 г. В частности, у западнокамчатской горбуши подходы уверенно растут в четные годы (рис. 2); у кеты р. Большой (одной из важнейших промысловых рек Западной Камчатки) промышленный вылов также значительно увеличивается в этот же период (За-варина, 2010); значительно возросли возвраты нерки (Дубынин, 2012). При этом снижения заполнения нерестилищ у этих видов не отмечается, что также свидетельствует в пользу второго фактора — недоучета производителей кижуча. Основная при-
1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014
Год
• Подход, N —— Пропуск на нерестилища, Б
Рис. 1. Динамика подходов кижуча и его пропуска на нерестилища на Западной Камчатке
200 180 160 140 S 120 х 100 2 80 60 40 20 0
1957 1960 1963 1966 1969 1972 1975 1978 1981 1984 1987 1990 1993 1996 1999 2002 2005 2008 2011
Год
-----Потомки, R ----------Пропуск на нерестилища, S
Рис. 2. Динамика подходов горбуши и ее пропуска на нерестилища на Западной Камчатке
I
\ А 1 i i А 1 Л 1
V\ Л А А /\ л А /\ 1 /
г Л А П 1\1\и /А /\ IaIa /\ /л\
\ Л А
■ У,,, Л У У Vy ,,,1, У, У, У,
чина такой ситуации заключается, конечно, в резком сокращении финансирования авиаучетных работ во второй половине 2000-х гг. Здесь приходится констатировать, что авиаучет кижуча—наиболее «больной» вопрос для мониторинга заходов лососей на Камчатке: сроки его нерестового хода самые поздние и растянутые. Как правило, в период проведения авиаучетных работ именно для кижуча нехватка финансирования имеет определяющее значение и обусловливает недоучет его производителей на нерестилищах. С 2006-2007 гг. общий уровень финансирования авиаучетных работ, и до этого времени дефицитный, снизился по отношению к предыдущему периоду в 1,5 раза, что резко сократило учет подходов кижуча в пользу более востребованных промыслом видов.
Таким образом, заниженная оценка числа зашедших в реки производителей влечет за собой достаточно важные последствия в оценке будущих
подходов в сторону занижения. Эту проблему можно решить двумя способами: первый — вносить экспертные корректировки в оценку родительского запаса, второй — использовать модели прогнозирования численности, не привязанные к этой величине.
Итак, цель данной работы — дать прогноз подхода лососевых, применяя несколько различных моделей, принадлежащих к разным классам. В качестве примера используются данные по популяции кижуча Западной Камчатки. Предполагается, что получив схожие либо одинаковые оценки будущих подходов, используя при этом несколько различных методов моделирования, мы будем иметь более обоснованную окончательную оценку подхода.
МАТЕРИАЛ И МЕТОДИКА
Исходным материалом для моделирования являлись следующие временные ряды (рис. 1-5):
80
10------------------------------------------------------------------------------------------------------------------------------------------------------
0 ----1---1---1---1--1---1--1---1---Г—I 1------1---1---1--1---1---1--1---1---1--1---1---1 I 1 I I I I I I I 1-------------------1—1-----1---1--1 I I
1970 1972 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008
Год
-----Доля возраста 2+
Рис. 3. Колебания долевого участия возрастной группы 2+ у кижуча
1957 1960 1963 1966 1969 1972 1975 1978 1981 1984 1987 1990 1993 1996 1999 2002 2005 2008 2011
Год
-----Потомков горбуши на ед. производителя R/S
Рис. 4. Колебания относительной величины R/S (количества потомков на производителя) у горбуши Западной Камчатки
- общего подхода кижуча (численность —
N Л
киж '
- пропущенных на нерест производителей кижуча (родители — $киж) и обусловленного ими пополнения (потомки — R );
4 киЖ'
- долевого участия возрастных когорт кижуча в пополнении (доли возрастов — С2+, С3+, С4+);
- соотношения родителей и потомков горбуши (коэффициент воспроизводства Rг0р^/Sг0рl).
Взаимосвязь исходных данных, методов разведочного анализа (автокорреляции и кросскорреляции), преобразования и регрессии, а также анализа ошибок, показана на схеме (рис. 6).
В работе рассматриваются четыре регрессионные модели прогнозирования подходов кижуча, которые можно разбить на два основных класса:
- модели, основанные на экстраполяции временного ряда,
- модели, основанные на теории зависимости пополнения от родительского нерестового запаса.
Из первого класса использовались модели авторегрессии и проинтегрированного скользящего среднего отклонений (ARIMA), с помощью которых строили модельный ряд подходов кижуча (Ы) и расщепления возрастных когорт кижуча в пополнении. Модель была предложена Боксом и
Год
Потомки кижуча, Я
Рис. 5. Пополнение кижуча Западной Камчатки в динамике
Рис. 6. Схема проведения исследований
Дженкинсом (1974) и включает в себя параметры авторегрессии и скользящего среднего. Имеется три типа параметров модели: параметры авторегрессии (р), порядок разности ^), параметры скользящего среднего (д). В обозначениях Бокса и Дженкинса модель записывается как АШМА (р, d, д). Например, модель (0, 1, 2) содержит 0 параметров авторегрессии (р) и 2 параметра скользящего среднего (д), которые вычисляются для ряда после взятия разности с лагом 1.
Для идентификации модели ARIMA (определения количества параметров) крайне важен предварительный анализ автокорреляций исходного и преобразованного с помощью взятия разности (дифференцированного) временного ряда. Автокорреляции представляют в виде автокорреляционной и частной автокорреляционной функций (АКФ и ЧАКФ).
АКФ показывает статистическую взаимосвязь между случайными величинами из одного ряда, но взятых со сдвигом по времени. Максимумы на автокорреляционной функции соответствуют циклическим колебаниям.
ЧАКФ является более углубленной версией обычной автокорреляционной функции. Ее отличительной особенностью является исключение корреляционной зависимости между наблюдениями внутри лагов, т. е. частная автокорреляционная функция на каждом лаге отличается от обычной автокорреляционной функции на величину удаленных автокорреляций с меньшими временными лагами. Следовательно, частная автокорреляционная функция более точно характеризует автокорреляционные зависимости внутри временного ряда.
Визуальный анализ АКФ и ЧАКФ помогает определить, какой процесс (авторегрессии или скользящего среднего) преобладает, а также количество параметров этих процессов. Если, например, АКФ имеет синусоидальный вид или экспоненциально убывает, то преобладает процесс авторегрессии. Если, наоборот, такой вид имеет ЧАКФ, то имеет место процесс скользящего среднего. В большинстве случаев имеют место смешанные процессы, а чтобы получить модель временного ряда, хватает одного-двух параметров авторегрессии и (или) одного-двух параметров скользящего среднего. Основные сочетания параметров р и д и вид соответствующих АКФ и ЧАКФ определены в (Бокс, Дженкинс, 1974).
Из класса моделей временных рядов, кроме ARIMA, рассматривается также модель распределенного лага (DLM), с помощью которой моделируется пополнение кижуча. Эта модель учитывает наличие разных промежутков времени между разными частями зависимой и независимой переменной, где в качестве ряда независимой переменной выступает ряд отношения между потомками и родителями горбуши №гор/$горе). Об основных предпосылках выбора этой объясняющей переменной для оценки пополнения кижуча будет сказано ниже, в подразделе, рассматривающем эту модель.
Из второго класса использовали достаточно популярные модели Рикера и Шепарда.
Аппроксимация модельных значений к фактическим проводилась путем минимизации суммы квадратов отклонений. После проводились оценка значимости уравнения регрессии (с помощью анализа дисперсий) и оценка значимости модельных коэффициентов (по /-критерию Стьюдента).
Кроме регрессионного анализа, для большинства моделей использовались промежуточные и разведывательные статистические методы, такие как: преобразование рядов, автокорреляционный и кросскорреляционный анализы.
Для удобства и ускорения расчетов использовались программы:
- MS Excel (предварительная аппроксимация, графика);
- Statistica 6.0 (регрессионный анализ моделей DLM, Рикера и Шепарда, анализ ошибок);
- EViews 8 (анализ автокорреляций, моделирование и оценка моделей ARIMA, графика).
Анализ отклонений модельных данных от фактических проводился по следующим показателям:
- имеется ли автокорреляция в отклонениях модельных данных от фактических: если автокорреляция наблюдается, значит какой-либо фактор не охватывается моделью и ее необходимо доработать;
- распределению отклонений модели от фактических данных по частотному ряду. Иными словами, анализируется вид распределения ошибок, проверяется гипотеза, принадлежит ли это распределение к нормальному.
Для подтверждения гипотезы о нормальном распределении ошибок использовались следующие критерии:
- Колмогорова-Смирнова с поправкой Лил-лиефорса ^ПИейэге, 1967), при расчете которого находится максимальное отклонение d между выборочной и теоретической интегральными функциями распределения и проверяется нулевая гипотеза о нормальном распределении выборки;
- Жарка-Бера (Тащие, Вега, 1981), в котором сравниваются коэффициенты асимметрии и эксцесса между выборочным и нормальным распределением. Чем ближе распределение ошибок к нормальному, тем меньше данная статистика отличается от нуля. При достаточно большом значении статистики вероятность нулевой гипотезы о нормальности распределения будет мала, и тогда будет основание ее отвергнуть.
Если в отклонениях модели отсутствует автокорреляция и частотное распределение ошибок близко к нормальному, то мы можем заключить, что модель в целом адекватно описывает фактический ряд наблюдений.
Изначально предполагалось, что те модели для одного и того же фактического ряда, которые будут признаны статистически значимыми и адекватными, должны дать достаточно близкие прогнозы; как минимум, их доверительные интервалы должны перекрываться.
Таким образом, следующим шагом анализа является необходимость построения доверительного интервала прогнозных значений на заданном уровне вероятности. Для этого рассчитывались стандартные ошибки точечных прогнозов, учитывающие как неопределенности самой модели регрессии (стандартная ошибка регрессии), так и неопределенности отдельных значений факторных переменных (Брюков, 2011):
зЕи=s ■ , (8)
где X — матрица исходных значений факторных переменных по всему временному ряду; -.
X — транспонированная матрица исходных значений факторных переменных по всему временному ряду;
Х{ — матрица-столбец значений факторных переменных для момента времени ^
— транспонированная матрица-столбец значений факторных переменных для момента времени t,
S — стандартное отклонение уравнения регрессии.
Доверительный интервал точечного прогноза будет равен:
а,_= У + !8ЕШ, (9)
где t — аргумент функции Лапласа, при а=0,05, t=1,96.
После вычисления доверительных интервалов прогнозных значений остается последний этап верификации модели — это сравнение плотности вероятности прогноза по ней с аналогичным показателем у других моделей. Если плотность вероятности прогноза такой модели не пересекается с другими плотностями (на заданном уровне доверия), то, следовательно, вес прогноза такой модели наименьший. В общем, мы проводим сравнительный анализ прогнозов по различным моделям и должны принять решение об окончательном прогнозе. Сделать это можно следующими способами:
- первый: найти средневзвешенный прогноз. т. е. сначала присвоить веса каждому прогнозу, затем найти их средневзвешенное значение;
- второй: перемножить плотности вероятности прогнозов и, следовательно, найти плотность вероятности того, что все прогнозы верны.
РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ
Авторегрессионная модель подходов кижуча (АЯ1МА). Идентификация модели авторегрессии и скользящего среднего начинается с анализа АКФ и ЧАКФ. Автокорреляционная функция исходного ряда имеет синусоидальный вид с полуперио-дом около 11 лет, частная автокорреляционная функция имеет значимую корреляцию на единичном лаге и не имеет на других (рис. 7).
Так как исходный временной ряд имеет восходящий тренд, преобразуем его в стационарный путем взятия разности (дифференцирования ряда) первого порядка Л1^=Nt-Nt_1, где Nt — ряд подходов кижуча в год t (обозначим как ряд D(N,1)).
Partial Correlation
АС РАС Q-Stat Prob
■I 1 1 ] 0,634 0,444 0,634 0,071 7,3149 11.181 0,007 0.004
1 ] 1 2
Ш ' 1 1 1 3 0,324 0.030 13.419 0,004
□ 1 1 с 1 4 0,131 -0,172 13.817 0,008
1 1 1 1 1 5 0,064 0,029 13,920 0,016
1 1 1 1 6 0,034 0,020 13,952 0,030
• с 1 • с 1 7 -0,142 -0,244 14,591 0,042
1 а 1 1 с 1 8 -0.290 -0,223 17,664 0,024
' СИ 1 1 1 9 -0,318 -0,019 21,960 0,009
' щ 1 1 1 10 -0.330 0,011 27,520 0,002
' CZ 1 1 с 1 11 0.353 -0,140 35,671 0.000
• с= 1 1 1 1 12 -0,301 -0,045 43,381 0,000
и получим также коррелограммы преобразованного ряда. Однако АКФ и ЧАКФ не имеют значимых корреляций (рис. 8).
Возьмем еще раз разность между соседними членами ряда D(N,1), получим ряд D(N,2XAW^=AW^-AW^_1) и проанализируем коррелограммы (рис. 9). На них уже присутствуют значимые корреляции на втором лаге. Фактически значимые результаты автокорреляций ряда разности второго порядка по сравнению аналогичными для ряда с разностью первого порядка говорят о наличии нелинейного тренда в исходном ряде подходов.
Так как вид АКФ и ЧАКФ до третьего лага мало отличаются друг от друга, а на следующих лагах корреляции уменьшаются и беспорядочны, то сложно сказать, какой процесс преобладает: авторегрессии второго порядка или скользящего среднего этого же порядка. Поэтому в дальнейшем анализе идентификация модели ЛКМЛ идет методом перебора моделей с количеством параметров авторегрессии р и скользящего среднего q, показанных в таблице 1.
В результате наиболее значимой моделью осталась модель авторегрессии второго порядка и скользящего среднего третьего порядка:
Л2 Nt = а х А2 Nt-2 + Ь х -3 + ,
где а — параметр авторегрессии второго порядка, Ь — параметр скользящего среднего ошибки третьего порядка, £— ошибка.
В EViews данная модель записывается как:
D(N,2) = AR( 2)+ МА(3).
Результаты регрессионного анализа и оценки значимости коэффициентов показаны в (табл. 2-4). Сравним фактические и модельные коррелограммы (рис. 10).
Таблица 1. Варианты моделей ЛКМЛ
tocorrclation Partial Correlation АС РАС Q-Stat Prob
1 1 1 • 1 0,061 0,061 0,0631 0,802
СГ 1 ' С= 1 2 -0,364 -0,369 2,5331 0,282
П 1 1 =1 I 3 0,190 0,281 3,2671 0.352
С 1 • 1 4 -0,084 -0,353 3,4269 0,489
1 1 С 1 5 -0,379 -0,154 7,0067 0,220
] ' 1 1 1 6 0,079 -0,040 7,1806 0,304
—1 1 • ZI 1 7 0,338 0,251 10,828 0,146
] 1 1 3 • 8 0,071 0,131 11,015 0.201
с 1 і с 1 9 -0,129 -0,129 11,759 0,227
1 • [ 1 10 -0,008 -0,083 11,762 0,301
с 1 ' с 1 1] -0,084 -0,134 12,295 0,342
с 1 1 ] < 12 -0,175 -0,076 15,731 0,204
Рис. 8. АКФ и ЧАКФ ряда первой разности подходов кижуча D(N,1)
Autocorrelation ' С
Partial Correlation • С
С
AC РАС Q-Stat Prob
і -0,139 -0,139 0,3162 0,574
2 -0,587 -0,618 6,4203 0,040
3 0,364 0,239 9,0003 0,029
4 0,210 -0,105 9,9520 0,041
5 -0,349 0,018 12,914 0,024
6 -0,129 -0,287 13,377 0,037
7 0,269 0,066 15,737 0,028
8 -0,005 -0,150 15,738 0,046
9 -0,162 0,165 17,010 0,049
10 0.121 -0,039 17,955 0.056
11 -0.084 -0,043 17,973 0,082
12 -0,175 -0,075 20,193 0.064
Рис. 9. АКФ и ЧАКФ ряда разности второго порядка подходов кижуча D(N,2)
0,8
■а 0,4
Я
і
О
0
о
3 -0,4 <
-0,8
с 0,4
о 7
.2
ё
о
о О
3 -0,4
Я
<2 0,8
/ _ 1 ■ • 1
1 4 / ГГ- і 1
1 13 4 5 6 7 8 9
/ • . " 1
\
\
1 2 3 4 5 ' 6 7 8 9
-----Actual -------Theoretical
Рис. 10. Фактические и модельные коррелограммы ряда разности второго порядка подходов кижуча D(N,2)
Таблица 2. Регрессионная статистика
0
1
2
+
+
+
+
Множ. R R-квадрат Норм. R-квадрат Стандартная ошибка Наблюдения
0,S36
0,697
0,665
0,1S5
11
3
+
+
+
Как видно из рисунков, теоретические корре-лограммы достаточно хорошо описывают фактические значения как АКФ, так и ЧАКФ.
Сравним визуально ряд и модельный ряд D (N,2), а также ряд ошибок на рис. 11. Модель достаточно хорошо описывает фактический ряд.
Далее суммируем (интегрируем) полученный модельный ряд D (N,2) два раза, это обратная операция к взятию разности ряда (дифференцированию). В итоге получаем модельный ряд подходов (рис. 12). Собственно, после этой процедуры модель ARMA превращается в модель ARIMA.
Анализ остатков показывает адекватность модели. Гипотеза о нормальном распределении остатков не отвергается, значимых автокорреляций в ряду остатков не наблюдается (рис. 13-14). В заключение рассчитаем стандартные ошибки индивидуальных прогнозов по формуле (8), с помощью которых определим доверительный интервал (95%) точечных прогнозов, в том числе и на 2014 год.
Запомним интервальный прогноз на 2014 год (2,59±0,43 млн экз.) и перейдем к следующей модели.
Модель распределенного лага (DLM). Применительно к задаче расчета пополнения методом распределенных лагов в общем виде модель можно запи-
жизнь других видов лососей и, в конечном виде, на их распределение и численность. Отношений «хищник-жертва» между различными видами тихоокеанских лососей не наблюдается, но так как горбуша служит основным источником биогенов в речной экосистеме, то влияние это не прямое, а косвенное, выражающееся через кормовую базу. Отнерестившиеся производители горбуши являются основным источником биогенов в речной экосистеме и через трофические связи оказывают Таблица 3. Дисперсионный анализ
df SS MS F Значимость F
Регрессия 2 0,755 0,377 10,925 0,004
Остаток 9 0,311 0,035
Итого 11 1,065
сать как:
Rt = A) x xt +А x Xt-1 + в2 x xt-2 + •••
---Residual-----Actual ---Fitted
Рис. 11. Фактический и модельный временные ряды разности второго порядка подходов кижуча D(N,2) и ошибки модели ARMA
(10,
+в х
где Rt—пополнение кижуча: количество потомков от производителей в год t, в0 п — коэффициенты,
Х1!п — независимые объясняющие переменные.
В качестве объясняющей переменной х был выбран коэффициент воспроизводства западнокамчатской горбуши R/S (отношение потомков к производителям, показывает величину потомков на одного производителя). Причины выбора такой переменной следующие.
1. Как численно доминирующий вид, горбуша оказывает большое влияние на
Таблица 4. Оценка параметров регрессии
Рис. 12. Фактический и модельный временные ряды (ЛЯІМЛ) подходов кижуча N с доверительным интервалом прогнозных значений (а=0,05)
Переменная Коэффициенты Стандартная ошибка t-статистика P-значение Нижние 95% Верхние 95%
AR(2) MA(3) -0,570 0,878 0,237 0,074 -2,399 11,798 0,040 8,9 X 10-7 -1,107 0,710 -0,032 1,047
влияние на выживаемость молоди тихоокеанских лососей на всех стадиях онтогенеза речного периода жизни (Бирман, 1985; Заварина, Шевляков, 2004; Шевляков, Заварина, 2004). Факты, свидетельствующие о влиянии горбуши на численность и прочие биологические показатели других видов лососей, отмечались неоднократно (Бирман, 1985; Бугаев, 1995, 2000; Bugayev, Dubynin, 2000;
Шунтов, Темных, 2004).
2. Чередование численности пополнения кижуча в четные и нечетные годы, наблюдающееся в последнее время (рис. 5, 2002-2008 гг.). Такой феномен можно объяснить влиянием смены у горбуши депрессивной линии нечетных лет и продуктивной линии четных лет.
3. Показатель R/S у горбуши можно представить как интегральный показатель, обобщающий условия окружающей среды. Условно его высокие значения соответствуют режиму благоприятных экологических факторов, а низкие — соответственно, неблагоприятных.
4. Продолжительность цикла жизни горбуши два года, что минимум на два года меньше, чем у кижуча (основные возрастные группы которого на Западной Камчатке 2+ и 3+). Это может дать заблаговременные оценки пополнения кижуча, а следовательно, даст возможность прогнозировать.
5. Кросскорреляционный анализ между временными рядами пополнения кижуча R и коэффициента воспроиз-
J киж ~ ~ 1
водства горбуши R/S (рис. 15) показывает значимые корреляции на лагах с запаздыванием и говорит о нетривиальных отношениях между двумя этими рядами.
Следующим шагом анализа является преобразование исходных данных по коэффициенту воспроизводства горбуши.
В некоторые годы (1992, 1999 и 2009) данный параметр (количество потомков на производителя) достигает неправдоподобно высоких значений (рис. 16).
Это можно объяснить недоучетом производителей на нерестилищах, особенно в малочисленные нечетные годы депрессивной линии горбуши. При не-
больших подходах производителей ошибка их учета нарастает экспоненциально. Подробно об этом говорится в работе (Шевляков, Маслов, 2011).
Рис. 13. Гистограмма остатков @=0,16953, Колмогорова-Смирнова тест: р>0,2; Лиллиефорс тест: р>0,2; Гипотеза о нормальности распределения остатков не отвергается. Вероятность нормального распределения — 47% по Жарку-Бера)
Рис. 14. АКФ и ЧАКФ остатков модели
Рис. 15. Кросскорреляционный анализ между временными рядами пополнения кижуча R и коэффициента воспроизводства горбуши R/S
Год
---1*/5 горбуши ----преобразованный ИУЯ
Рис. 16. Исходный и преобразованный ряды коэффициента воспроизводства горбуши R/S
Было принято решение преобразовать ряд Я/5 горбуши с помощью экспоненциального сглаживания, по формуле:
г =
Г + г,-1 к + г, _21 к 3
(11)
где г — коэффициент воспроизводства горбуши
ад
г — преобразованный коэффициент воспроизводства горбуши,
k — параметр, определяющий вес переменной г на лагах -1 и -2. На нулевом лаге — к0 = 1. Таким образом, со временем его значение изменяется в геометрической прогрессии. Этот параметр определяется одновременно с коэффициентами Ьп при минимизации суммы квадратов отклонений приведенного ниже уравнения (12).
Итоговая модель зависимости пополнения кижуча от коэффициента воспроизводства горбуши: Rtкиж = Ь0 х г + Ь хГ- + Ь2 х г_2 +... (12)
+Ь х г
П t-П
Преобразованный ряд коэффициентов воспроизводства горбуши при оптимальном для итоговой модели к показан на рис. 16. При минимизации квадратов отклонений (ошибок) модели от реальных данных на различных временных отрезках, выяснилось, что модель, начиная с 1993 года, достаточно хорошо описывает временной ряд численности поколений кижуча (рис. 17).
На высоком уровне значимости оказались коэффициенты на лагах -12, -10, -9, -7 и 0. Остальные коэффициенты незначимы, и соответствующие лаги были исключены из модели. Регрессионная статистика и значимость уравнения регрессии показаны в таблицах 5-7. Данные показатели очень высокие. Причина несоответствия модели фактическому ряду до 1993 г. заключается в следующем: в 1978 г. был запрещен японский дрифтерный промысел в 200-мильной экономической зоне (Кур-мазов, 2001), вследствие чего, как минимум до этого года, данные по уловам
Таблица 7. Оценка параметров регрессии
горбуши и кижуча (следовательно, и данные по подходам к местам нереста) являются неполными. Однако если из 1993 г. вычесть максимальный лаг объясняющей переменной (-12), мы получим 1981 год. Вспомним, что преобразование ряда коэффициента воспроизводства горбуши шло методом сглаживания по трем последним годам. Иными словами, объясняющая переменная для 1981 года г1981 несет в себе информацию и за предыдущие два года. Вычитая их, мы получим 1979 год, первый год с новыми условиями промысла и отличающейся статистикой.
На рис. 18 показаны фактические и смоделированные данные общей численности пополнения кижуча Западной Камчатки с 95%-м доверительным интервалом. Численность пополнения кижуча от родителей 2010 г. составит К2010 = 1,64 млн экз., а от родителей 2011 г. — Я2011 = 2,395 млн экз.
Моделирование процентных долей в возрастных когортах потомков кижуча. Моделирование численности потомков — это только половина
Таблица 5. Регрессионная статистика
Множ. R R-квадрат Норм. R-квадрат Стандартная ошибка Наблюдения
0,998
0,992
0,893
0,072
15
Таблица 6. Дисперсионный анализ
ат SS MS F Значимость F
Регрессия 6 10,637 1,773 311,355 6,5х10-10
Остаток 9 0,051 0,006
Итого 15 10,689
Год
----Я факт. Я модел.
Рис. 17. Фактическая и модельная численность поколений кижуча Западной Камчатки (модель DLM)
Переменная Коэффициенты Стандартная ошибка ^статистика Р-значение Нижние 95% Верхние 95%
-12 0,109 0,019 5,783 0,000 0,066 0,152
-10 0,117 0,015 7,954 0,000 0,084 0,150
-9 0,058 0,008 7,002 0,000 0,039 0,077
-7 0,047 0,009 5,346 0,000 0,027 0,067
0 -0,028 0,006 -4,638 0,001 -0,042 -0,014
к 1,351 0,202 6,688 0,000 0,894 1,808
решения задачи прогнозирования численности будущего подхода рыбы. Необходимо также определить доли возрастных групп в будущем подходе, и сделать это можно различными путями. Одним из самых простых путей является взятие средних значений долей за последние несколько лет. Такой подход является оправданным, когда с течением времени доли возрастных групп мало изменчивы и остаются на стабильном уровне. В противном случае, когда за рассматриваемый период мы наблюдаем достаточно сильные колебания возрастных долей, а также можем предположить цикличность, лучше воспользоваться моделированием возрастного состава. У кижуча Западной Камчатки нерестовое стадо преимущественно (в среднем 98% от подхода за период с 1970 по 2009 года) состоит из особей возрастов 2+ и 3+ (соответственно третий и четвертый года жизни) (Зорбиди, 2010), поэтому в нашем случае достаточно смоделировать временной ряд какой-либо одной группы. Возьмем ряд изменения доли возраста 2+
(рис. 3). После визуального анализа можно предположить, что ряд представляет собой не случайные колебания. Действительно, мы видим разные режимы: двухлетние циклы в 1974-1978, 1990-1997 гг.. которые чередуются с 5-, 6-летними циклами в 1979-1989 гг. и 3-, 4-летними циклами в 1997-2009 гг. Поэтому для моделирования возьмем ряд, начиная с 1997 г. Используем модель АШМА.
Так же, как и в показанной выше модели подходов АШМА, сначала приведем ряд к стационарному, путем взятия разности первого порядка и анализа АКФ и ЧАКФ (рис. 19).
Автокорреляции имеют значимые уровни на втором лаге, и их вид весьма напоминает аналогичные автокорреляции для ряда подходов кижуча, что неслучайно, ведь доля возрастной группы 2+ (как и других возрастных групп) — это доля именно от его подхода. В данном случае нельзя определенно сказать, какой процесс преобладает: авторегрессии 2-го порядка или скользящего сред-
Таблица 10. Оценка параметров регрессии
него 2-го порядка, поэтому так же, как и при моделировании подходов, для идентификации модели необходимо рассмотреть варианты моделей, показанные в таблице 1.
В итоге наиболее подходящей адекватной моделью, отвечающей условиям стационарности и обратимости, остается модель вида: D(С2+,1)=AR(2)+MA(1).
Статистическая значимость уравнения регрессии и коэффициентов, а также анализ ошибок, показаны в таблицах 8-10.
Таблица 8. Регрессионная статистика Множественный Я 0,9027
Я-квадрат 0,8148
Стандартная ошибка 7,08
Наблюдения___________________10__________________
Таблица 9. Дисперсионный анализ I df I SS I MS I F
I Значимость F
Регрессия 2 1339,84 669,92 13,356 0,003
Остаток 8 401,28 50,16
Итого 10 1741,12___________________________
Рис. 18. Фактическая и модельная численность поколений кижуча Западной Камчатки (модель DLM) с доверительным интервалом прогнозных значений (а=0,05)
Рис. 19. АКФ и ЧАКФ ряда первой разности доли возрастной группы 2+ ряд D(С2+,1)
Переменная Коэффициенты AR(2) -0,8537
МА(1)_________-0,9371
Стандартная ошибка 0,2280 0,0577
1-статистика | Р-значение Нижние 95%Верхние 95% -3,745 0,0057 -1,3794 -0,3280
-16,233 0,0000 -1,0703 -0,8040
Сравним фактические и модельные коррелограммы, а также фактический и модельный временные ряды долевого участия возрастной группы 2+ (рис. 20-21).
Анализ ошибок этой модели (рис. 2223) показывает возможность опровержения гипотезы о нормальном распределении остатков на уровне а=0,1 (по Лил-лиефорсу). Другие тесты эту гипотезу не опровергают. Автокорреляционный анализ тоже «на грани»: довольно высокое значение частной функции на 4-м лаге, хотя и не выходящее за пределы доверительного интервала.
Теперь, используя модель DLM пополнения, выполненную на первом этапе анализа, и модель ARIMA процентных долей в расщеплении поколений на возраста, выполненных на втором этапе, можно спрогнозировать общий подход кижуча Западной Камчатки в 2014 г. Для этого надо перемножить смоделированные значения пополнения и соответствующие им доли возрастов:
N =С2+ X R +С3+ X R
2014 2011 2011 2010 2010
Фактический подход, его модельные значения и 95%-й доверительный интервал индивидуальных прогнозов покажем на графике (рис. 24).
Анализ ошибок не опровергает гипотезу о нормальном распределении остатков. Вероятность гипотезы о нормальном распределении по тесту Жарка-Бера р=0,997 (рис. 25-26).
Как и в предыдущей модели, запомним полученный результат: прогноз подхода кижуча в 2014 году — 2,18±0,27 млн экз.
Следующие две модели, которые мы рассмотрим, принадлежат к классу моделей зависимости пополнения от нерестового запаса.
Модель Рикера. В качестве уравнения, описывающего зависимость пополнения от запаса, используется формула Рикера (2).
Как и в предыдущей модели DLM, численность подхода рассчитывается с помощью модельных значений пополне-
1,0
с
В 0,5
я
Е ■о о =
<
-1.0
\ _ \
1 \ /' N /
1 2 3 4 5 6 7
— Actual Theoretical
Рис. 20. Фактические и модельные коррелограммы ряда D(C2+,1)
Рис. 21. Фактический и модельный ряды изменения доли возрастной группы 2+ у кижуча
-10 -5
Рис. 22. Гистограмма остатков @=0,24656, Колмогорова-Смирнова тест: р>0,2; Лиллиефорс тест: р<0,1; гипотеза о нормальности р пределения отвергается на уровне а=0,1; по тесту Жарка-Бе гипотеза о нормальности распределения не отвергается)
ас-
ера
tocorrelaiion Partial Correlation АС РАС Q-Stat Prob
С 1 1 С 1 1 -0,176 -0,176 0,4146
1 1 1 1 1 2 -0,418 -0.464 3.0403
Z! 1 1 1 1 3 0,230 -0,052 3,9438 0,047
С 1 1 1 1 4 -0,241 0,483 5,1039 0,078
с 1 1 С 1 5 -0,121 -0,222 5,4571 0.141
=1 1 1 С 1 6 0,334 -0,146 8,8121 0,066
В 1 1 □ 1 7 0,131 0,193 9,4958 0,091
с 1 1 1 1 8 -0,174 -0.041 11,304 0,079
I 1 1 1 9 -0,064 0,008 11,800 0,107
ния R, возрастов 2+ и 3+ и их долей в подходе (используется модель АЛ1МА возрастного состава):
4 = ^\-зХ Л-з + С+_4х ^. (13)
Основной трудностью при моделировании стал оптимальный выбор интервала временных рядов нерестового запаса и пополнения кижуча. Разброс точек в облаке данных на диаграмме зависимости Л от S (рис. 27) достаточно большой для выбора однозначной и статистически значимой модели. Даже для промежутка данных в последние десять лет (2000-2009 гг.) параметры модели будут явно не стационарными (рис. 28). Поэтому набор исходных данных для регрессионного анализа пришлось значительно сократить до последних пяти лет наблюдений (2004-2009 гг.). Обусловлено это скачкообразным изменением параметров модели, начиная с 2008 года. На рис. 28 хорошо видно, что линии регрессии № 1-3 изменяются незначительно, однако следующие две линии № 4 и № 5, характерные включением данных 2008 и 2009 годов соответственно, значительно отличаются от предыдущих. Максимум куполообразной кривой за два года вырастает почти в два раза, иными словами, почти вдвое повышается эффективность воспроизводства.
Результаты регрессионного анализа 4
показаны в таблицах 11-13. Несмотря на
„ „ з
короткий временной ряд исходных данных величиной В ПЯТЬ лет, модель И ее 2 параметры а и Ь получились статистически значимыми с высокой степенью достоверности.
Модельные значения подхода на фоне фактических данных показаны на рис. 29. Напомним, что фактические значения 2004-2008 гг. не учитывались в модели, а модельные значения для этих лет показаны с помощью обратной экстраполяции. Для этого промежутка видим, что модельные значения всегда выше фактических, однако стоит заметить, что тенденции возрастания и убывания численности практически не изменяются.
Ограниченный пул исходных данных существенно расширяет границы доверительного интервала. Однако с другой
Таблица 11. Регрессионная статистика
Множ. R R-квадрат Стандартная ошибка Наблюдения 0,899 0,809 0,225 5
Таблица 12. Дисперсионный анализ
SS df MS F Значимость F
Регрессия 8,5109 2 4,2555 83,992 0,0023
Остаток 0,1519 3 0,0506
Итого 8,6629 5
Mean -0,006220
Median -0,002069
Maximum 0,111964
Minimum -0,129517
Sld.Dcv. 0,060159
Skewness -0,019508
Kurtosis 3,078806
Год
-----N набл. N молел. —........С195%.
Рис. 24. Фактический и модельный ряды подходов кижуча (модель DLM), с доверительным интервалом прогнозных значений (а=0,05)
Scries: SHR1ES01 Sample 1 15 Observations 15
Jarque-Bera
Probability
0,004833
0,997586
-0,15 -0,10 -0,05 0.0 0,05 0,10
Рис. 25. Г истограмма остатков @=0,15664, Колмогорова-Смирнова тест: р>0,2; Лиллиефорс тест: р>0,2; гипотеза о нормальности распределения остатков не отвергается. Вероятность нормального распределения — 99% по Жарку-Бера)
Autocorrelation □
al Correlation АС РАС Q-Stat Prob
□ 1 1 0,203 0,203 0,7505 0,386
С 1 2 -0,097 -0,144 0,9346 0.627
С 1 3 -0,246 -0,207 2.2238 0,527
=) 1 4 0,111 0.216 2.5113 0.643
1 5 0.132 0,020 2,9531 0,707
С 1 6 -0,117 -0,217 3,3379 0.765
1 7 -0,414 -0,304 8,7883 0.268
• 8 -0,169 -0,013 9,8298 0,277
с 1 9 -0,010 -0,125 9,8344 0,364
1= 1 10 -0,014 -0,210 9,8440 0,454
с 1 И -0,174 -0,132 11,777 0,381
п 1 12 0.044 0.169 11,942 0.450
стороны, анализ ошибок не отвергает гипотезу о нормальности распределения остатков, и модель можно признать адекватной (рис. 30-31).
С этими оговорками примем к дальнейшему рассмотрению прогноз подхода кижуча на 2014 год, составивший 1,72±0,62 млн экз.
Модель Шепарда. Следующей моделью из класса «запас-пополнение» рассмотрим модель Шепарда, вида:
R=aS/(b+S2). (14)
В отличие от модели (5), мы предлагаем принять параметр с постоянно равным 2 (константе), исходя из следующих соображений. Во-первых, мы понизим количество параметров в модели, а тем самым повысим ее статистическую значимость, ведь наблюдений у нас всего пять. Во-вторых, такую модель (как и модель Рикера) можно будет представить в линейном виде, что облегчает регрессионный анализ и вычисление необходимых статистик. И в-третьих, есть возражения теоретического плана: если параметр с не равен целым 0, 1 или 2 (например, с=3,65), то невозможно аналитически объяснить, — что означает нерестовый запас S в степени 3,65. Однако рассуждения на эту тему выходят за рамки данной работы.
А численность будем рассчитывать, как и в предыдущей модели, по формуле (13). Результаты регрессионного анализа и анализа ошибок представлены в таблицах 14-16 и на рис. 32-33. Статистическая значимость уравнения регрессии выше, чем у модели Рикера, однако параметр Ь статистически незначим. С другой стороны, доверительный интервал индивидуальных прогнозных значений несколько уже, чем по модели Рикера. Как и в предыдущем случае, при анализе отклонений модели, гипотезы о нормальном распределении остатков не отвергаются.
Таблица 13. Оценка параметров регрессии
2
1,8
1,6
1,4
1,2
і
| 0.8 с 0,6 0,4 0,2 0
• •
••
О
0,2
1
1,2
0,4 0,6 0.8
Родители 8, млн экз.
Рис. 27. Пополнение кижуча в зависимости от пропущенных на не рест рыб (1971-2009 гг.)
Рис. 28. Данные по родителям и потомкам кижуча Западной Камчатки за период 2000-2009 гг. и проведенные к ним кривые Рикера за пятилетние промежутки
----N факт.------N модел. -----С195%
Рис. 29. Фактический и модельный ряд подходов кижуча (по модели Рикера), с доверительным интервалом прогнозных значений (а=0,05)
Переменная Коэффициенты а 17,468
Ь 3,814
Стандартная ошибка 3,916 0,552
^статистика Р-значение Нижние 95% Верхние 95% 4,461 0,021 5,006 29,929
6,906 0,006 2,056 5,572
Таблица 14. Регрессионная статистика
Таблица 15. Дисперсионный анализ
Множ. R 0,931 ат SS MS F Значимость F
R-квадрат 0,866 Регрессия 8,556 2 4,278 120,582 0,0014
Стандартная ошибка 0,188 Остаток 0,106 3 0,036
Наблюдения 5 Итого 8,663 5
Таблица 16. Оценка параметров регрессии
Переменная Коэффициенты а 0,6080
Ь 0,0228
Стандартная ошибка 0,1137 0,0199
1-статистика | Р-значение |Нижние 95%|Верхние 95% 5,347 0,0128 0,2461 0,9698
1,146 0,3348 -0,0405 0,0862
Рис. 30. Гистограмма остатков @=0,2235, Колмогорова-Смирнова тест: р>0,2; Лиллиефорс тест: р>0,2; гипотеза о нормальности распределения остатков не отвергается — 80% по Жарку-Бера)
Рис. 31. АКФ и ЧАКФ остатков модели 3 ■
Рис. 32. Гистограмма остатков @=0,2551, Колмогорова-Смирнова тест: р>0,2; Лиллиефорс тест: р>0,2; гипотеза о нормальности распределения остатков не отвергается — 69% по Жарку-Бера)
Рис. 33. АКФ и ЧАКФ остатков модели
Модельные значения подходов кижуча в сравнении с фактическими данными показаны на рис. 34. Прогноз на 2014 г. по этой модели составляет 2,16±0,52 млн экз.
РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Получив прогнозы по четырем рассмотренным моделям, определим, какая модель является лучшей. Для этого сравним их статистики в сводной таблице 17.
Практически по всем статистикам, лучшей является модель БЬМ. С прогнозом по этой модели с точностью до 0,2 млн экз. совпадает прогноз, данный по модели Шепарда, однако она имеет более широкий интервальный прогноз и худшие статистики, что, как мы помним, обусловлено небольшим набором наблюдений, однако по статистикам она в основном лучше, чем ее «одноклассник», модель Рикера. Наихудший результат по стандартным ошибкам регрессии и прогноза показывает модель Рикера, кроме того, она единственная показывает снижение численности подхода в 2014 г., по сравнению с 2013 г.
Чтобы дать средневзвешенную оценку интервального прогноза на 2014 год, необходимо найти средневзвешенные оценки точечного прогноза и его стандартного отклонения. Для этого будем ранжировать модели по 4-балльной системе, где наивысший балл будет получать модель с наилучшей статистикой. Балльные оценки показаны в табл. 18.
В результате наивысший балл набирает модель БЬМ (28 баллов), затем следуют модели Шепарда (16 баллов) и ARIMA (15 баллов), наименьший результат у модели Рикера (11 баллов). Набранное количество баллов будет служить весовым коэффициентом вклада каждой модели. Рассчитаем средневзвешенный прогноз:
N
ш . N + Ж . N + Ж
гг ЛШМЛ у ЛШЫЛ^ гг DLЛ у DLЛ^ ггРикер
. N + Ж N
Рикер Шепард Шепард
срвзв.
(15)
Ж +Ж +Ж +Ж
гг ЛШМЛ^ гг DLЛ^ ' гРикер ^ ' Шепард
где N — прогнозы, а Ж — весовые коэффициенты для каждой модели. Получаем:
15 • 2,59 + 28 • 2,18 +11-1,72 +16 • 2,16
N
срвзв.
15 + 28 + 11 +16
Далее вычисляем средневзвешенное стандартное отклонение:
срвзв.
1
У (Ж - N )2
/ у ' і Срвзв.1
= 2,19млн экз.
Затем средневзвешенную стандартную ошибку (где п — количество моделей):
, 0,2642
БЕ =
= 0,1321
і =1
п -1
где п — общее количество баллов,
4,8162 70 -1
= 0,2642
Рис. 34. Фактический и модельный ряд подходов кижуча (по модели Шепарда), с доверительным интервалом прогнозных значений (а=0,05)
4п л/4
Находим доверительный интервал при а=0,05: С10=0,13211,96=0,259 Следовательно, средневзвешенный прогноз подхода кижуча на 2014 год составит:
= 2,19 ± 0,259 млн экз.
2014 7 7
Другим путем обобщения полученных результатов является представление полученных интервальных прогнозов в виде плотностей вероятности. Имея точечные прогнозы и их доверительные интервалы С10 95, предполагаем, что точечный прогноз представляет собой среднее, приближенное к математическому ожиданию нормального распределения его ошибок, и представляет собой первый параметр нормального распределения. Вторым параметром является стандартное отклонение от среднего, и в нашем случае оно равно стандартной ошибке прогноза (табл. 17).
Таблица 17. Основные статистики рассмотренных моделей
Модель Статистическая значимость уравнения регрессии Станд. ошибка регрессии Я2 Станд. ошибка прогноза Макс. р Н0 коэф. Информ. критерий Шварца ВІС Информ. критерий Акаике АІС
F Р
ЛЯІМЛ 10,92524 0,0040 0,1858 0,70 0,1896 0,0310 -3,1306 -4,0137
БЬМ 311,355 6,5Х10-10 0,0716 0,99 0,1407 0,0012 -4,5959 -30,619
Рикер 83,9920 0,0023 0,2251 0,81 0,3177 0,0210 -2,8496 0,7227
Шепард 120,5815 0,0014 0,1884 0,87 0,2658 0,3348 -3,2058 -1,0586
Таблица 18. Набранные моделями баллы
Модель Стат. значимость уравнения регрессии Станд. ошибка регрессии Я2 Станд. ошибка прогноза Стат. значимость коэф. регрессии ВІС ЛІС Всего
ЛЯІМЛ 1 3 1 3 2 2 3 15
БЬМ 4 4 4 4 4 4 4 28
Рикер 2 1 2 1 3 1 1 11
Шепард 3 2 3 2 1 3 2 16
Используя эти два параметра, мы можем получить функции нормального распределения и представить их как плотности вероятности прогнозов по каждой модели (рис. 35).
Учитывая, что все полученные прогнозы являются только различными (отчасти совпадающими) оценками одного и того же реального будуще-
го подхода, можно выделить область, где все четыре оценки будут учитываться; иными словами, все данные прогнозы сбудутся. В теории вероятность одновременного наступления событий определяется как произведение вероятностей каждого из них. Графически на рис. 36 произведение вероятностей прогнозов показано вместе с плотностью
0,35
0,30
0,25
0,20
0,15
0,10
0,05
0,50 0,65 0,80 0,95 1,10 1,25 1,40 1,55 1,70 1,85 2,00 2,15 2,30 2,45 2,60 2,75 2,90 3,05 3,20 3,35
N. млн экз.
| АШМА ■ РИКЕР ШЕПАРД ■ ОЬА Рис. 35. Плотности вероятности прогнозов подхода кижуча Западной Камчатки на 2014 г.
!
■
■МІРІ 11
1.7 1,8 1,9 2,0 2,1 2,2 2,3 2,4 2,5 2,6 2,7
N. млн экз.
■ произведение вероятностей средневзвешенный прогноз
Рис. 36. Обобщенные прогнозы по четырем моделям: как произведение плотностей вероятности и как средневзвешенный прогноз
вероятности ранее вычисленного средневзвешенного прогноза. А прогноз, определенный как произведение вероятностей четырех моделей, составит (при а=0,05):
М.т=2,21±0,125 млн экз.
2014
В итоге можно резюмировать, что при применении двух способов обобщения результатов по рассмотренным моделям мы получили приблизительно один и тот же точечный прогноз (с различием в 0,02 млн экз.) подхода кижуча в 2014 году: 2,18-2,20 млн экз. При этом доверительные интервалы отличаются: средневзвешенный доверительный интервал практически в два раза шире, чем аналогичный у произведения плотностей вероятности. Объясняется это тем, что произведение плотностей вероят-
ности имеет более жесткие ограничительные рамки (все прогнозы верны и сбудутся). Средневзвешенный же плотностной прогноз, напротив, допускает возможность того, что один из двух прогнозов по крайним моделям неверен (модели АШМА или Рикера) и является более подходящим при осторожной политике управления рыбными запасами.
ЗАКЛЮЧЕНИЕ
Предложенный подход к решению задачи прогнозирования численности тихоокеанских лососей не является строго регламентированным. В нем должны соблюдаться основные моменты: применение нескольких прогностических моделей, обязательный анализ отклонений модельных данных от фактического материала, сравнение моделей между собой и выбор окончательного прогноза. Во всем остальном исследователю даются широкие возможности по применению различных видов прогностических моделей, методов разведывательного анализа, процедур сглаживания и преобразования временных рядов, выбора статистических критериев и методов аппроксимации, определения выбросов из данных или длины временного ряда. Во всех этих случаях не обойтись также без должного опыта и определенной интуиции. С другой стороны, мы против применения в прогнозе какой-либо одной модели без сравнения ее результатов с показателями альтернативных моделей. Также интересно будет узнать, насколько полученный нами прогноз подхода западнокамчатского кижуча оправдается.
СПИСОК ЛИТЕРАТУРЫ
Бабаян В.К. 1988. Математические методы теории рыболовства (модели изолированных популяций) // Рыб. хоз-во. Обзоры по информационному обеспечению общесоюзных научно-технических программ. М.: ЦНИИТЭИРХ. 76 с.
Бабаян В.К. 2000. Предосторожный подход к оценке общего допустимого улова (ОДУ): Анализ и рекомендации по применению. М.: ВНИРО, 192 с. Баранов Ф.И. 1918. К вопросу о биологических основаниях рыбного хозяйства // Изв. отдела рыбоводства и науч.-промысл. Исследований. № 1 (1). С. 84-128.
Бивертон Р., Холт С. 1969. Динамика численности промысловых рыб. Пер. с англ. М.: Пищ. пром-сть. 248 с.
Бирман И.Б. 1985. Морской период жизни и вопросы динамики стад тихоокеанских лососей. М.: Агропромиздат. 208 с.
Бокс Дж., Дженкинс Г. 1974. Анализ временных рядов, прогноз и управление. М.: Мир. Вып. 1-2. 604 с.
Брюков В.Г. 2011. Как предсказать курс доллара. Эффективные методы прогнозирования с использованием Excel и EViews. М.: ЦИПСиР. 272 с. Бугаев В.Ф. 1995. Азиатская нерка (пресноводный период жизни, структура локальных стад, динамика численности). М.: Колос. 464 с.
Бугаев В.Ф. 2000. Мир реки. Петропавловск-Камчатский: Северная Пацифика. 32 с.
Дубинин В.А. 2012. Об оптимуме производителей нерки на нерестилищах бассейна р. Озерной в современный период / Матер. Всерос. науч. конф., посвящ. 80-летнему юбилею ФГУП «КамчатНИРО». Петропавловск-Камчатский: КамчатНИРО. С. 302-308.
Заварина Л.О., Шевляков Е.А. 2004. Возможный механизм формирования цикличности урожайных поколений кеты на северо-восточном побережье Камчатки / Сохранение биоразнообразия Камчатки и прилегающих морей. Матер. V науч. конф. Петропавловск-Камчатский. С. 52-55.
Заварина Л.О. 2010. О динамике биологических показателей и тенденциях изменения численности кеты (Oncorhynchus keta) р. Большой (ЮгоЗападная Камчатка) // Исслед. водн. биол. ресурсов Камчатки и сев.-зап. части Тихого океана: Сб. науч. тр. КамчатНИИ рыб. хоз-ва и океанографии. Вып. 18. С. 38-57.
Запорожец О.М., Шевляков Е.А., Запорожец Г.В., Антонов Н.П. 2007. Возможности использования данных о нелегальном вылове тихоокеанских лососей для реальной оценки их запасов // Вопр. рыболовства. Т. 8. № 3 (31). С. 471-483.
Зорбиди Ж.Х. 2010. Кижуч азиатских стад. Петропавловск-Камчатский: КамчатНИРО. 306 с. Криксунов Е.А. 1995. Теория пополнения и интерпретация динамики популяций рыб // Вопр. ихтиологии. Т. 35. № 3. С. 302-321.
Курмазов А.А. 2001. Международно-правовые условия освоения морских биологических ресурсов в Тихом океане // Мир. океан: использование биол. ресурсов. М.: ВИНИТИ. С. 24-41.
Хилборн Р., Уолтерс К. 2001. Количественные методы оценки рыбных запасов. Выбор, динамика, неопределенность. СПб.: Политехника. 228 с.
Шевляков E.A., Заварина Л.О. 2004. Об особенностях динамики численности и методиках прогнозирования запасов кеты Oncorhynchus keta W. (Salmonidae) Западной Камчатки // Исслед. водн. биол. ресурсов Камчатки и сев.-зап. части Тихого океана: Сб. науч. тр. КамчатНИИ рыб. хоз-ва и океанографии. Вып. 7. С. 181-186.
Шевляков E.A., Маслов A.B. 2011. Реки, определяющие воспроизводство тихоокеанских лососей на Камчатке, как реперы для оценки заполнения нерестового фонда // Изв. Тихоокеан. НИИ рыб. хоз-ва и океанографии. Т. 164. С. 114-139. Шевляков E.A., Дубинин B.A., Зорбиди Ж.Х., Заварина Л.О., Попова T.A., ApтюхинаН.Б., Горин С.Л., Коваль О.О. 2013. Современное состояние лососевого реки Большой (Западная Камчатка): воспроизводство, промысел, управление // Изв. Тихоокеан. НИИ рыб. хоз-ва и океанографии. Т. 174. С. 3-37. Шунтов В.П., Темных О.С. 2004. Взгляд на лососевую путину 2004 через призму итогов изучения и промысла лососей в 2003 г. // Рыбн. хоз-во № 2. С.26-27.
Anderson E.D. (Ed.). 2002. The Raymond J.H. Be-verton lectures at Woods Hole, Massachusetts. Three Lectures on Fisheries Science Given May 2-3, 1994. U.S. Dep. Commer., NOAA Tech. Memo. NMFS-F/ SPO-54. 161 p.
Bugayev V.F., Dubynin V.A. 2000. Factors influencing abundance of sockeye salmon (Oncorhynchus nerka) from the Ozernaya River, Southwest Kamchatka // Recent Changes in Ocean Production of Pacific Salmon J.H. Helle, Y. Ishida, D. Noakes and V. Radchenko (ed.). North Pac. Anadromous Fish Com. Bull. 2. Vancouver. Canada. P. 181-189.
Cushing D.H. 1973. The dependence of recruitment on parent stock // J. Fish. Res. Board. Canada. V. 30. P. 1965-1975.
Deriso R.B. 1980. Harvesting strategies and parametr estimation for an age. Structured model // Can. J. Fish. Aquat. Sci. V. 37. P. 268-282.
Hoff J.C. 1983. A practical guide to Box-Jenkins forecasting. Belmont, Calif.: Lifetime Learning Publications. 316 p.
Jarque C.M., Bera A.K. 1981. Efficient tests for normality, homoscedasticity and serial independence of regression residuals: Monte Carlo evidence // Economics Letters. V. 7 (4). P. 313-318.
Lilliefors H. 1967. On the Kolmogorov-Smirnov test for normality with mean and variance unknown // Journal of the American Statistical Association/ V. 62. № 318 (Jun., 1967). P. 399-402.
Pankratz A. 1983. Forecasting with univariate Box-Jenkins models: Concepts and cases. N.Y.: Wiley. 562 p. Paulik G.J. 1973. Studies of the possible form of the stock-recruitment curve // Rapp. P. - V. Reun. Cons. Perm. Intern. Explor. Mer. V. 164. P. 302-315.
Quinn T.J. II. 2003. Ruminations on the development and future of population dynamics models in fisheries // Natural Resource Modeling 16 (4). P. 341-392. Ricker W. 1954. Stock and recruitment // Journal of the fisheries research board of Canada. V. 11. № 5. P. 559-623.
Sheperd J.G. 1982. A versatile new stock-recruitment relationship for fisheries, and the construction of sustainable yield curves // J. Cons. Intern. Explor. Mer. V. 40 (1). P. 67-75.
Vandaele W. 1983. Applied time series and Box-Jenkins models. London: Academic Press. 417 p.
