Desarrollo del indicador de máxima superficie en riesgo de incendio

La variable aleatoria “superficie anual afectada por incendios forestales” (y) se define como la suma de la superficie afectada por todos los incendios ocurridos a lo largo de un año en una unidad territorial y, como toda variable aleatoria, posee una distribución de densidad f(y) y una función de distribución F(y). Dicha variable se conoce en la literatura anglosajona como “área quemada anual” (annual burned area) y no debe confundirse con la variable “superficie afectada por incendio” (SAI) en una unidad territorial, la cual es la más usada comúnmente para modelar distribuciones de tamaño de incendios. Las distribuciones de frecuencia de la variable SAI exhiben asimetría positiva, similar a la que se ilustra en la Figura 1 para la variable aleatoria y; sin embargo, se caracterizan por tener una cola muy pesada (^{Jiang, Zhuang, Flannigan, & Little, 2009}; ^{Katz, Brush, & Parlange, 2005}), que dificulta el modelaje, incluso con el uso de distribuciones de valores extremos (^{Alvarado, Sandberg, & Pickford, 1998}; ^{Jiang & Zhuang, 2011}).

Figura 1 Variable aleatoria y (superficie anual afectada por incendios forestales) con distribución de densidad f(y). Se ilustra la máxima superficie en riesgo de incendio (SeR) con probabilidad máxima (P _max ) de 5 %. 

En lo que resta del documento, la variable aleatoria y será referida como la superficie afectada y equivale a la suma de superficies afectadas por incendios a lo largo de un año en una unidad territorial. Cabe señalar que esta variable es el indicador más usado para calificar la severidad de una temporada de incendios, el riesgo de la superficie bajo aprovechamiento maderable o para asociar la calidad del aire y salud en presencia de incendios (^{Wiitala & Carlton,1994}). La variable también tiene una asimetría positiva (Figura 1), pero el peso de la cola derecha es menor, debido a que los incendios individuales con afectaciones muy grandes son un sumando más en la suma de superficies afectadas en un año y no representan un valor único en el extremo de la distribución, como sucede en distribuciones de eventos individuales.

Si f(y) es conocida, entonces la máxima superficie en riesgo de incendio (SeR) en una unidad territorial, dado un nivel de probabilidad P _max , se define como la superficie que corresponde al valor de y en la cola derecha de f(y) para un nivel de probabilidad P _max ; es decir, SeR = F ^-1 (1- P _max ). La SeR es un indicador análogo al valor en riesgo (^{Jorion, 2000}) conocido popularmente como VaR y usado como medida de riesgo en el sector financiero. La SeR se ilustra en la Figura 1 con P _max = 5 %. Este indicador, medido en superficie, tiene el atributo de aislar los eventos de valores muy altos (superficie total afectada en un año) de los que ocurren en 1- P _max de las ocasiones. Con ello se logra estimar la magnitud de una pérdida poco frecuente, pero de gran escala, sin el efecto de los eventos (superficie afectada en un año) más comunes.

La determinación del valor de P _max depende del riesgo que se esté dispuesto a tomar en la estimación de la superficie de máxima afectación. El VaR financiero usa un nivel de probabilidad de 5 %, aunque ocasionalmente utiliza 2.5 % y 1 % (^{Dowd, 2003}). Para la estimación de la SeR, el valor P _max = 5 % parece apropiado, ya que aproxima la frecuencia relativa de incendios de gran magnitud (^{Taylor et al., 2013}) y también puede asociarse a un periodo de retorno (PR = 11-Fy ) de 20 años.

Otros indicadores relacionados

Los incendios forestales de gran magnitud regularmente se asocian al cambio climático, mayor acumulación de combustibles o al manejo deficiente de incidentes (^{Fernandes, Pacheco, Almeida, & Claro, 2016}); por ello, los incendios siguen un patrón no estacionario derivado de la misma característica de sus principales causales, lo cual origina que su distribución de probabilidades cambie de periodo a periodo (^{Read & Vogel, 2016}). Este cambio en la distribución de incendios grandes y, por consecuencia, de las superficies afectadas entre el año t-k y el año t, se puede evaluar por la diferencia en el valor de la SeR en ambos tiempos. A esta medida se le denominará en lo subsecuente SeR incremental (ISeR); su cálculo para un periodo de k años, dada una P _max (ISeRk,Pmax) se realiza como ISeRk,Pmax=SeRt, Pmax-SeRt-k, Pmax , donde SeR _i,j corresponde a la SeR en el periodo i calculada con una P _max = j.

La ISeR brinda una métrica de la magnitud y el sentido en que el riesgo de ocurrencia de una superficie afectada cambia en el periodo de k años. Un aumento en la ISeR en un periodo de k años indica que el riesgo aumenta a tasa creciente (no lineal), dada la no linealidad de la cola de la distribución de frecuencias (^{Dowd, 2003}).

La razón entre la superficie afectada por incendios en un año y la superficie forestal total, en cada unidad territorial, es una expresión relativa de la superficie afectada respecto a la superficie forestal. Esta proporción constituye una nueva variable aleatoria x (x = superficie total incendiada/superficie forestal), la cual también tiene una función de densidad f(x). De manera análoga a la SeR, la “máxima proporción de superficie forestal con riesgo de ser afectada por incendios” (PSeR) con un nivel de confianza P _max , se definió como la proporción en la cola derecha de f(x) dado un nivel de probabilidad P _max .

Si f(x) es conocida para dos momentos consecutivos, es posible calcular una diferencia en la proporción máxima de superficie en riesgo de incendio (PSeR) en un periodo de k años y dada una P _max (ISeRk,Pmax). Esta diferencia se denominará PSeR incremental (IPSeR), la cual se estima como IPSeRk,Pmax=PSeRt,  Pmax-PSeRt-k, Pmax , donde PSeR _i,j corresponde a la PSeR en el periodo i calculada con una P _max = j.

Datos

La superficie afectada por incendios forestales por entidad y por año se obtuvo de la CONAFOR para el periodo 1970-2018. La superficie con cobertura forestal varía año con año derivado de los procesos de cambio de uso del suelo. La estimación de esta cobertura tiene un nivel variable de precisión. Por ello, la superficie forestal se consideró constante en el tiempo y se determinó con base en la última cobertura de vegetación disponible (serie VI del ^{Instituto Nacional de Estadística y Geografía [INEGI], 2015}). Este supuesto hace que la forma de la distribución empírica de la variable PSeR [f(x)] sea similar a aquella de la variable SeR [f(x)], dado que X = Y/k, donde k es constante en todos los años.

Métodos

Cada una de las series de datos de superficie afectada a causa de incendios forestales, por entidad, se ajustó a varios modelos de función de distribución para identificar el modelo de mejor ajuste. Los modelos probados fueron: gama, exponencial, lognormal y Weibull, selección realizada con base en las recomendaciones de varios autores (^{Alvarado et al., 1998}; ^{Cui & Perera, 2008}; ^{Cumming, 2001}; ^{Jiang & Zhuang, 2011}). La función de probabilidad se ajustó a partir de la estimación de una distribución empírica [ F~y ] (^{Forbes, Evans, Hastings, & Peacock, 2011}). Esta distribución se ajustó al modelo de prueba (gama, exponencial, lognormal y Weibull) con el uso del estadístico Dmax de la prueba Kolmogorov-Smirnov, como criterio de minimización de desviaciones. El estadístico se define como Dmax=supyF~y-F^y , donde F^y representa la función de distribución predicha bajo el modelo hipotético. Esta metodología aseguró que los estimadores cumplieran con el criterio Dmax mínimo desde el ajuste. La estimación se realizó con la aplicación Risk Simulator™, desarrollada para Excel.

El desempeño de cada modelo se evaluó con la significancia del estadístico Dmax y dos criterios cualitativos. El primer criterio calificó con un valor de 1 al modelo con mejor significancia, 2 al modelo con la siguiente mejor significancia y así sucesivamente hasta 4 que corresponde al modelo de menor significancia. La calificación de cada modelo en todas las entidades se sumó y se usó como un primer indicador del modelo más recomendable. El segundo criterio calificó con 1 al modelo que brindó el mejor ajuste y 0 para los demás. La suma de estos valores, dividido entre el número de entidades, define la proporción de ocasiones en que un modelo brindó el mejor ajuste.

Una vez elegido el modelo de mayor significancia se procedió a calcular la SeR para valores de P _max de 1 %, 5 %, 20 % y 90 % para las series 1970-2008, 1970-2013 y 1970-2018, siguiendo el modelo de mejor ajuste para cada entidad (i. e. no se usó un modelo genérico para todas las entidades). Estos valores de P _max se eligieron dado que corresponden a periodos de retorno de 100, 20, 5 y 1 año, respectivamente.

A fin de identificar el poder predictivo de los indicadores desarrollados, la predicción para el año t+1 de un modelo compuesto por la SeR y la ISeR se contrastó con la predicción para el mismo año basada en la variable superficie afectada en el año anterior (SAA), variable que en teoría tiene mayor poder predictivo. El modelo de predicción con la SeR y la ISeR se basa en un modelo lineal de dos componentes: i) un estimador de riesgo sistemático o inevitable relacionado con el riesgo inherente de la unidad territorial y ii) un estimador de la tendencia de la superficie en riesgo, el cual como se ha señalado, aproxima la acumulación de combustibles. El componente de riesgo sistemático se definió por la SeR con P _max = 90 % (SeR_0.90), debido a que esta superficie es la esperada en un periodo de retorno de un año, mientras que el componente de tendencia se aproximó con la ISeR de los últimos k = 5 años (ISeR_5,0.90). El modelo de predicción tiene la forma básica SAA _t+1 = β0 + β1 SeRt,  0.90 + β2 ISeR5, 0.90 + e, donde las ß _i corresponden a parámetros del modelo y e corresponde al término de error. Este modelo es análogo a un algoritmo de búsqueda en el cual, dado un punto de referencia, se aproxima el siguiente punto con una dirección y un tamaño de paso. El tamaño de paso está definido por el valor de la SeR _Pmax , mientras que la dirección se aproxima por la ISeR _k,Pmax estimada con la misma P _max . Variantes de este modelo se compararon con el modelo de predicción de un retraso, definido por SAA _t+1 = β0 + β1 SAA _t + e. La comparación se realizó con el criterio de R².

Forma funcional de f(y)

Los ajustes tuvieron significancia alta para el estadístico Dmax en todos los estados, a excepción de Baja California Sur, entidad en la que fue necesario eliminar las observaciones con nula superficie afectada para mejorar la calidad de los ajustes. El modelo lognormal brindó la mejor calificación de acuerdo con los criterios cualitativos, seguido del modelo gama, exponencial y al final el modelo Weibull (Cuadro 1). Lo anterior muestra que la distribución de incendios sigue una función lognormal en la mayoría de los casos, tal y como lo han referido varios autores (^{Holmes et al., 2008}; ^{Reed, 2001}; ^{Strauss, Bednar, & Mees, 1989}).

Superficie máxima en riesgo de incendio (SeR)

La Figura 2 muestra la superficie máxima en riesgo de incendio (SeR) para valores de P _max al 1 % (SeR_{1 %}) y al 5 % (SeR_{5 %}). En este gráfico sobresalen los estados de mayor superficie forestal (Chihuahua, Durango, Jalisco, Quintana Roo, Oaxaca, Chiapas y Guerrero), derivado de su mayor exposición al riesgo (mayor superficie arbolada). Algunos estados como Coahuila, Sonora, Baja California, Zacatecas y Nuevo León mostraron una SeR significativamente grande; sin embargo, el resultado puede ser engañoso dado que la superficie reportada como afectada no necesariamente está relacionada con la superficie arbolada.

Figura 2 Estimaciones de superficie máxima en riesgo de incendio (SeR) para valores de probabilidad máxima (P _max ) al 1 y 5 %. 

La Figura 2 resalta las diferencias entre la SeR_{1 %} y la SeR_{5 %} para estados como Quintana Roo, Coahuila, Sonora, Oaxaca, Campeche, Chiapas y Querétaro. El diferencial muestra una mayor ocurrencia de grandes afectaciones en periodos de retorno de 100 años (SeR_{1 %}). Presumiblemente, el diferencial podría estar asociado a condiciones particulares de acumulación de combustibles, condiciones climáticas (temperaturas extremas y vientos específicos), regímenes de fuego, eficiencia de las acciones de supresión o incluso de condiciones socioeconómicas que crean un entorno de mayor riesgo para la propagación de eventos de gran escala. Estas diferencias contrastan con entidades como la Ciudad de México, Estado de México o Sinaloa, donde es evidente que los eventos extremos (i. e. afectaciones totales grandes en un año) en un periodo de retorno de 100 años, no son proporcionalmente más grandes en relación con aquellos que podrían observarse en un periodo de 20 años.

Máxima proporción de la superficie en riesgo de incendio (PSeR)

A diferencia de la SeR, la PSeR es comparable entre unidades territoriales, independientemente de su tamaño y cantidad de superficie forestal; por ello es un índice relativo. La PSeR evaluada con P _max = 5 % y P _max = 1 % en cada entidad, para el periodo 1970-2018, se muestra en la Figura 3. En esta figura, considerando P _max = 5 %, los valores de las entidades con PSeR menor de 3 % se dividieron de aquellas con mayor proporción. Los resultados mostraron algunos patrones interesantes. Los estados de mayor superficie forestal se encuentran dentro del grupo de entidades con PSeR < 3 %; esto es, las afectaciones proporcionalmente mayores no ocurren en las entidades de mayor superficie forestal. De esta forma, entidades como Coahuila, Baja California, Nuevo León y Tabasco, con baja proporción de superficie forestal, presentan altos valores de PSeR. Asimismo, el caso de Baja California destaca con valores muy altos de PSeR (26 %). No obstante, este valor es engañoso dado que la PSeR es relativa a la cobertura forestal y gran parte de los incendios en esta entidad corresponden a pastizales o vegetación semiárida. En una situación similar se encuentran los estados de Coahuila y Nuevo León.

Figura 3 Estimaciones de la máxima proporción de la superficie en riesgo de incendio (PSeR) para valores de probabilidad máxima (P _max ) de 1 % y 5 %. Entidades con PSeR < 3 % considerando P _max = 5 % (a) y PSeR ≥ 3 % considerando P _max = 5 % (b), en el periodo 1970-2018. 

Otro aspecto relevante en los resultados es la PSeR con P _max = 1 % (PSeR_{1 %}) para Quintana Roo, entidad con alta incidencia de eventos grandes; su PSeR_{1 %} es proporcionalmente más alta que la PSeR_{5 %}, lo que sugiere mayor probabilidad de que se observen eventos más aislados de gran magnitud que eventos de mediana magnitud. El mismo caso se observa para entidades como Campeche, Tamaulipas, Zacatecas, Yucatán, Querétaro y Baja California.

Proporción de la superficie en riesgo incremental (IPSeR)

Para ejemplificar la PSeR incremental (IPSeR) se usó un periodo k = 10 años. La diferencia se obtuvo comparando la PSeR del periodo 1970-2018 con aquella del periodo 1970-2008 (i. e. el periodo evaluado es 2009-2018). En todos los casos, la comparación se hizo con la forma funcional f(x) que tuvo mejor ajuste para cada estado. La Figura 4 muestra tres grupos de entidades: el primero compuesto por entidades donde la IPSeR es creciente, el segundo grupo con entidades donde prácticamente no hay cambio, y el tercer grupo donde la IPSeR es decreciente para el periodo 2009-2018. Los aumentos en la IPSeR superiores a 2 % solo se observan en la Ciudad de México, Baja California, Quintana Roo y Aguascalientes. Por el contrario, reducciones en la IPSeR solo se observan en Coahuila; el resto de las entidades parece tener cambios muy pequeños durante el periodo analizado.

Figura 4 Valores de proporción de la superficie en riesgo incremental (IPSeR_{10, 5 %}) de incendios entre 2009 y 2018. 

La IPSeR es un indicador de la dirección en que la PSeR cambia durante el periodo analizado y muy probablemente se relaciona con la cantidad de combustible disponible o la mayor incidencia de agentes causales en una unidad territorial. Un ejemplo de su relación con la disponibilidad de combustibles son los estados de Coahuila, Chiapas, México y Nuevo León, que tuvieron superficies incendiadas importantes en el periodo previo (2009-2018) y causaron una reducción en el volumen de combustibles. Tal reducción presumiblemente disminuyó las afectaciones posteriores a este periodo, lo cual se refleja en una IPSeR < 0. Por el contrario, entidades con IPSeR positivos (Ciudad de México, Baja California, Quintana Roo y Aguascalientes) se caracterizan por haber presentado superficies extensas afectadas en un periodo mayor de 10 años, antes de 2009, por lo que muy probablemente recuperaron el inventario de combustibles en el periodo evaluado (2009-2018) y muestran mayor riesgo de tener superficies extensas afectadas después de dicho año. Un resultado que refuerza este argumento es que, en estas entidades, el número promedio de incendios en el periodo analizado fue estadísticamente superior al promedio de toda la serie.

Predicción de superficies afectadas

El Cuadro 2 muestra los resultados de los ajustes realizados para la predicción de la SAA en el año 2014, a partir de variantes de los modelos señalados en la sección metodológica. La ISeR corresponde a aquella entre los periodos 2009-2013 (k = 5). En todos los casos se usó una transformación logarítmica, debido a que el modelo de mejor ajuste es un modelo lognormal. El modelo C brindó estimadores cuya suma de ß _i es cercana a la unidad, lo cual es consistente con el modelo teórico SAA _t+1 = β0 + β1 SeRt,  0.90 + β2 ISeR5, 0.90 + e y una ß ₀ restringida a cero. Por el contrario, en el modelo A no se restringió ß ₀ = 0, aunque el ajuste mostró que este estimador es estadísticamente igual a cero (ß = 0), lo que confirmó el resultado esperado para el modelo 3. El modelo D incluyó un componente no lineal para la SeR_0.90 (doble logaritmo de SeR_0.90) y muestra que la suma de ß _i se alejó de la unidad, pero con estimadores más robustos. Finalmente, el modelo autorregresivo (B) mostró un ajuste extraordinario para la predicción de 2014. No obstante, este modelo tiene una evidente restricción para predecir eventos poco frecuentes, ya que requiere de una calibración año con año, a diferencia de los modelos basados en SeR que son más consistentes. En resumen, el modelo D mostró el mejor ajuste; sin embargo, el modelo C presentó mejor parsimonia y fiabilidad.