La magnitud de la sequía

Cuando tienes una serie temporal al uso (sin entrar a definir qué es eso), uno puede aplicar descomposiciones tmabién al uso, como stl, para extraer tendencia y estacionalidad, de la forma

como en esta entrada previa.

Lluvia.

La serie de la lluvia es otra cosa. Uno ve si llueve o no llueve (típicamente no). Lo que uno no ve es la probabilidad, que varía a lo largo del año, de que llueva. Pasa lo mismo con monedas (uno ve caras o cruces, no la probabilidad de cara), clientes que compran (uno ve si compra o no, no la probabilidad de compra), etc. Pero la probabilidad existe y en estimarla consiste frecuentemente el trabajo de algunos.

En el caso de la lluvia, el que llueva un día $t$ (más de 1 mm, que es equivalente a 1 l/m²), o X_t en lo que sigue, es X_t \sim \text{Bernoulli}(p_t) y se puede suponer que las p_t tienen estacionalidad anual (en abril, aguas mil) y una tendencia histórica. Como arriba, pero con una diferencia: las p_t no se observan diariamente; a lo más, se observa si llueve o deja de llover un día determinado.

El modelo que propongo es

X_{m} \sim \text{Binom}(l(m), p_m)

donde X_{m} es el número de días que llueve (más de 1 mm) en el mes m; l(m) es el número de días del mes y p_m es la probabilidad de lluvia en cada uno de los días del mes. Esta probabilidad es

p_m = \frac{\exp(\eta_m)}{1 + \exp(\eta_m)}

donde \eta_m se descompone como la suma de una componente estacional (12 valores) y una tendencia global. Que son, según mis cálculos,

y

Se aprecia la sequía de 2005 y la de 2011 y 2012. Y la actual, por supuesto.

(Me asalta una duda: ¿contendrá esa tendencia final tan apocalíptica algún artefacto de mi opción de modelado?)

La probabilidad de lluvia (diaria) ha evolucionado así:

Los datos de pluviosidad histórica están en AEMET, pero sacarlos de ahí es tarea imposible. La AEMET está gobernada por funcionarios que dicen ñí. Así que los he bajado de NOAA; los de NOAA son funcionarios también, solo que estadounidenses y con vocación de servicio público. Como yo también la tengo, los comparto.

También el código, para referencia de todos. La magia, en todo caso, es producto del maravilloso paquete INLA de R. La parte más relevante del código es el lugar donde defino el modelo. La más discutible, es donde extraigo las estimaciones de los parámetros (me quedo con la moda de la posteriori). Pero es mejorable con poco esfuerzo.

Comenta

Your email address will not be published.

Puedes usar estas etiquetas y atributos de HTML: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code class="" title="" data-url=""> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre class="" title="" data-url=""> <span class="" title="" data-url="">