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 $latex X_t$ en lo que sigue, es $latex X_t \sim \text{Bernoulli}(p_t)$ y se puede suponer que las $latex p_t$ tienen estacionalidad anual (en abril, aguas mil) y una tendencia histórica. Como arriba, pero con una diferencia: las $latex 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 $latex X_{m}$ es el número de días que llueve (más de 1 mm) en el mes $latex m$; $latex l(m)$ es el número de días del mes y $latex 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 $latex \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.

library(ggplot2)
library(lubridate)
library(INLA)
library(plyr)

lluvia <- read.csv("1140486.csv")
lluvia <- lluvia[, c("DATE", "PRCP")]
lluvia$DATE <- as.Date(lluvia$DATE)

plot(lluvia$DATE, lluvia$PRCP, type = "l", xlab = "",
     ylab = "precipitación",
     main = "Precipitación diaria en Madrid\n(2000-2017)")

lluvia$llueve <- 1 * (lluvia$PRCP > 0)
lluvia$fecha <- as.numeric(lluvia$DATE)
lluvia$fecha <- 0 + lluvia$fecha - min(lluvia$fecha)
lluvia$doy   <- 1 + lluvia$fecha %% 365

# mensual

mensual <- lluvia
mensual$mes <- month(mensual$DATE)
mensual$ano <- year(mensual$DATE)

mensual <- ddply(mensual, .(ano, mes), summarize,
                    dias_lluvia = sum(llueve),
                    dias_totales = length(mes))

mensual$t <- 1:nrow(mensual)

formula <- dias_lluvia ~ f(mes, model="rw2", cyclic = TRUE, param = c(1,0.0001)) +
     f(t, model = "rw2", cyclic = FALSE, param = c(1,0.0001))
modelo <- inla(formula, family = "binomial", Ntrials = dias_totales, data = mensual, verbose = TRUE)

periodicidad <- sapply(modelo$marginals.random$mes, function(x) x[which.max(x[,2]),1])
tendencia <- sapply(modelo$marginals.random$t, function(x) x[which.max(x[,2]),1])
tmp <- modelo$marginals.fixed[[1]]
indep <- tmp[which.max(tmp[,2]),1]

prob <- indep + tendencia + rep(periodicidad, 100)[1:length(tendencia)]
prob <- exp(prob) / (1 + exp(prob))


png("/tmp/estacionalidad_lluvia_madrid.png", height = 400, width = 800)
plot(1:12, periodicidad, type = "l", xlab = "mes del año", ylab = "",
     main = "componente estacional de la probabilidad de lluvia\nMadrid, 2000-2017")
dev.off()

png("/tmp/tendencia_lluvia_madrid.png", height = 400, width = 800)
fechas <- seq(min(lluvia$DATE), max(lluvia$DATE), by = "1 month")
plot(fechas, tendencia, type = "l", xlab = "fecha", ylab = "",
     main = "tendencia de la probabilidad de lluvia\nMadrid, 2000-2017")
dev.off()

png("/tmp/probabilidad_lluvia_madrid.png", height = 400, width = 800)
plot(fechas, prob, type = "l", xlab = "fecha", ylab = "",
     main = "probabilidad diaria de lluvia\nMadrid, 2000-2017")
dev.off()