3 Formulación del Modelo Matemático
3.1 Modelo SDDE en tiempo continuo
3.1.1 Motivación y contexto teórico
Las ecuaciones diferenciales estocásticas con retraso (SDDE, por sus siglas en inglés) constituyen una extensión natural de las ecuaciones diferenciales estocásticas (SDE) clásicas, incorporando explícitamente la dependencia del pasado en la dinámica del sistema (Mao, 2007). Este tipo de modelos es particularmente relevante en sistemas biológicos donde los procesos fisiológicos presentan memoria temporal inherente, como es el caso de la fenología vegetal (Mohammed, 1984).
La necesidad de incorporar términos de retraso en modelos matemáticos de sistemas dinámicos fue reconocida inicialmente por Volterra en el contexto de ecuaciones integro-diferenciales (Volterra, 1931). Sin embargo, fue con el desarrollo de la teoría de procesos estocásticos y el Cálculo de Itô que se estableció un marco riguroso para el análisis de SDDEs (Itô, 1951).
3.1.2 Definición formal del modelo
Consideremos un espacio de probabilidad filtrado \((\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t \geq 0}, \mathbb{P})\) que satisface las condiciones usuales de completitud y continuidad derecha (Revuz & Yor, 1999). Sea \(W(t)\) un proceso de Wiener estándar (movimiento browniano) adaptado a la filtración \(\{\mathcal{F}_t\}\).
El modelo propuesto para la dinámica del índice de verdor GCC se define mediante la siguiente ecuación diferencial estocástica con retraso:
\[ dX(t) = \left[ \alpha + \beta X(t) + \gamma X(t-\tau) + \boldsymbol{\delta}^T \mathbf{C}(t) \right] dt + \sigma X(t) \, dW(t), \quad t \geq 0 \tag{3.1} \]
donde las variables y parámetros se definen en el contexto biológico de la siguiente manera:
\(X(t) \in \mathbb{R}^+\): Proceso estocástico del GCC, que representa la biomasa verde activa en el tiempo \(t\).
\(\alpha \in \mathbb{R}\): Tasa de crecimiento basal, asociada al crecimiento intrínseco del sistema en ausencia de forzantes externos.
\(\beta \in \mathbb{R}\): Coeficiente de retroalimentación actual, el cual modela el efecto de autorregulación del crecimiento.
\(\gamma \in \mathbb{R}\): Coeficiente de retroalimentación retardada, que representa la memoria fisiológica del sistema, vinculada a la translocación de fotoasimilados.
\(\tau \in \mathbb{R}^+\): Tiempo de retraso, correspondiente al periodo característico de los procesos fisiológicos involucrados.
\(\mathbf{C}(t) \in \mathbb{R}^3\): Vector de forzantes climáticas, definido como \(\mathbf{C}(t) = [\text{Rad}_t, \text{Temp}_t, \text{Prec}_t]^T\).
\(\boldsymbol{\delta} \in \mathbb{R}^3\): Vector de sensibilidades climáticas, que cuantifica la respuesta del sistema ante variaciones en las variables climáticas.
\(\sigma \in \mathbb{R}^+\): Intensidad del ruido estocástico, que representa la variabilidad ambiental no explicada por el modelo determinista.
\(W(t) \in \mathbb{R}\): Proceso de Wiener estándar, utilizado para modelar el ruido blanco gaussiano que perturba el sistema.
3.1.3 Condiciones iniciales y espacio de fases
A diferencia de las SDEs clásicas, las SDDEs requieren una función inicial definida en un intervalo pasado, debido a la dependencia del término \(X(t-\tau)\) (Mohammed, 1984). El problema de valor inicial se especifica como:
\[ X(t) = \phi(t), \quad \forall t \in [-\tau, 0] \tag{3.2} \]
donde \(\phi: [-\tau, 0] \to \mathbb{R}^+\) es una función continua y determinista que representa el historial del sistema antes del tiempo inicial \(t=0\).
El espacio de fases para una SDDE es infinito-dimensional, específicamente el espacio de Banach \(C([-\tau, 0], \mathbb{R})\) de funciones continuas del intervalo \([-\tau, 0]\) en \(\mathbb{R}\), equipado con la norma del supremo:
\[ \|\phi\| = \sup_{\theta \in [-\tau, 0]} |\phi(\theta)| \tag{3.3} \]
Esta característica fundamental distingue a las SDDEs de las SDEs y tiene implicaciones profundas para el análisis de existencia, unicidad y estabilidad de soluciones (Mao, 2007).
3.1.4 Propiedades matemáticas fundamentales
3.1.4.1 Existencia y unicidad de soluciones
Para garantizar la existencia y unicidad de soluciones a la SDDE (3.1), imponemos las siguientes condiciones estándar sobre los coeficientes (Mao, 2007):
Condición de Lipschitz global: Existe una constante \(K > 0\) tal que para todo \(x_1, x_2, y_1, y_2 \in \mathbb{R}\) y \(t \geq 0\):
\[ |\mu(t, x_1, y_1) - \mu(t, x_2, y_2)| + |\sigma(t, x_1) - \sigma(t, x_2)| \leq K(|x_1 - x_2| + |y_1 - y_2|) \tag{3.4} \]
donde \(\mu(t, x, y) = \alpha + \beta x + \gamma y + \boldsymbol{\delta}^T \mathbf{C}(t)\) es la función de deriva y \(\sigma(t, x) = \sigma x\) es la función de difusión.
Condición de crecimiento lineal: Existe una constante \(C > 0\) tal que para todo \(x, y \in \mathbb{R}\) y \(t \geq 0\):
\[ |\mu(t, x, y)|^2 + |\sigma(t, x)|^2 \leq C(1 + |x|^2 + |y|^2) \tag{3.5} \]
Teorema 3.1 (Existencia y unicidad). Bajo las condiciones (3.4) y (3.5), y dado un historial inicial \(\phi \in C([-\tau, 0], \mathbb{R}^+)\), existe una única solución global \(X(t)\) para todo \(t \geq 0\) que satisface la SDDE (3.1) y la condición inicial (3.2).
Demostración: Ver (Mao, 2007), Teorema 6.2, o Anexo D para una demostración detallada adaptada a nuestro modelo específico.
3.1.4.2 Verificación para el modelo propuesto
Para nuestro modelo específico, verificamos las condiciones anteriores:
Lipschitz: La función \(\mu(t, x, y) = \alpha + \beta x + \gamma y + \boldsymbol{\delta}^T \mathbf{C}(t)\) es lineal en \(x\) e \(y\), por lo tanto satisface la condición de Lipschitz con constante \(K = \max(|\beta|, |\gamma|)\).
Crecimiento lineal: Dado que \(\mu\) es lineal y \(\sigma(t, x) = \sigma x\) también es lineal, ambas satisfacen trivialmente la condición de crecimiento lineal.
Por lo tanto, el modelo propuesto garantiza existencia y unicidad de soluciones para cualquier historial inicial continuo \(\phi\).
3.1.5 Interpretación biológica de los términos
3.1.5.1 Término de deriva determinista
El término de deriva en la ecuación (3.1) puede descomponerse en cuatro componentes:
\[ \mu(t, X(t), X(t-\tau)) = \underbrace{\alpha}_{\begin{matrix} \text{Crecimiento} \\ \text{basal} \end{matrix}} + \underbrace{\beta X(t)}_{\begin{matrix}\text{Retroalimentación}\\ \text{actual}\end{matrix}} + \underbrace{\gamma X(t-\tau)}_{\begin{matrix}\text{Memoria} \\ \text{fisiológica}\end{matrix}} + \underbrace{\boldsymbol{\delta}^T \mathbf{C}(t)}_{\begin{matrix}\text{Forzantes} \\ \text{climáticas}\end{matrix}} \tag{3.6} \]
Crecimiento basal (\(\alpha\)): Representa la tasa de crecimiento intrínseca del sistema en ausencia de retroalimentación y forzantes externos. En el contexto vegetal, esto corresponde al crecimiento mínimo sostenido por procesos metabólicos básicos (Taiz & Zeiger, 2010).
Retroalimentación actual (\(\beta X(t)\)): Modela la autorregulación del crecimiento basada en el estado actual del sistema. Un valor \(\beta > 0\) indica crecimiento autocatalítico (aceleración), mientras que \(\beta < 0\) representa saturación o competencia intraespecífica (Diepenbrock, 2000).
Memoria fisiológica (\(\gamma X(t-\tau)\)): Captura explícitamente la dependencia del pasado, representando procesos como la translocación de fotoasimilados, acumulación de reservas o respuesta hormonal retardada (Thomas & Ougham, 2014). El signo de \(\gamma\) determina si el efecto es positivo (persistencia) o negativo (agotamiento de recursos).
Forzantes climáticas (\(\boldsymbol{\delta}^T \mathbf{C}(t)\)): Incorpora la influencia directa de variables ambientales sobre el crecimiento. Cada componente \(\delta_i\) cuantifica la sensibilidad del sistema a la variable climática correspondiente (Monteith & Moss, 1977).
3.1.5.2 Término de difusión estocástica
El término de difusión \(\sigma X(t) \, dW(t)\) modela la variabilidad ambiental no explicada por las variables climáticas incluidas explícitamente. La forma multiplicativa \(\sigma X(t)\) implica que la intensidad del ruido es proporcional al tamaño del sistema, lo cual es biológicamente razonable: sistemas más grandes (mayor GCC) exhiben mayor variabilidad absoluta (Kloeden & Platen, 1992).
3.1.6 Solución formal y representación integral
La solución de la SDDE (3.1) puede expresarse en forma integral utilizando el Cálculo de Itô (Øksendal, 2003):
\[ X(t) = X(0) + \int_0^t \left[ \alpha + \beta X(s) + \gamma X(s-\tau) + \boldsymbol{\delta}^T \mathbf{C}(s) \right] ds + \sigma \int_0^t X(s) \, dW(s) \tag{3.7}\label{eq-integral_solution} \]
Esta representación integral es fundamental para el análisis teórico y la discretización numérica del modelo. La integral estocástica \(\int_0^t X(s) \, dW(s)\) se interpreta en el sentido de Itô, lo cual implica que el integrando \(X(s)\) es evaluado en el extremo izquierdo de cada subintervalo de partición (Itô, 1951).
3.1.7 Momentos y propiedades estadísticas
Aunque la solución exacta de la SDDE (3.1) no puede expresarse en forma cerrada debido al término de retraso, podemos derivar ecuaciones para sus momentos estadísticos. Tomando esperanza en ambos lados de \(\ref{eq-integral_solution}\) y utilizando que la esperanza de la integral de Itô es cero:
\[ \mathbb{E}[X(t)] = X(0) + \int_0^t \left[ \alpha + \beta \mathbb{E}[X(s)] + \gamma \mathbb{E}[X(s-\tau)] + \boldsymbol{\delta}^T \mathbb{E}[\mathbf{C}(s)] \right] ds \tag{3.8} \]
Diferenciando respecto a \(t\) obtenemos una ecuación diferencial funcional determinista para la media:
\[ \frac{d}{dt} \mathbb{E}[X(t)] = \alpha + \beta \mathbb{E}[X(t)] + \gamma \mathbb{E}[X(t-\tau)] + \boldsymbol{\delta}^T \mathbb{E}[\mathbf{C}(t)] \tag{3.9} \]
Esta ecuación revela que la dinámica de la media del proceso estocástico sigue una DDE determinista con la misma estructura de retraso que el modelo original (Mohammed, 1984).
3.1.8 Consideraciones sobre la positividad del modelo
Dado que \(X(t)\) representa el índice de verdor GCC, que físicamente debe ser no negativo (\(X(t) \geq 0\)), es crucial verificar que el modelo preserve esta propiedad. Para SDEs con coeficientes lineales como la nuestra, se puede demostrar que si la condición inicial satisface \(\phi(t) \geq 0\) para todo \(t \in [-\tau, 0]\), entonces \(X(t) \geq 0\) casi seguramente para todo \(t \geq 0\) (Mao, 2007).
Esta propiedad de positividad es esencial para la interpretación biológica del modelo y garantiza que las simulaciones numéricas produzcan resultados físicamente realistas.
3.2 Discretización Numérica y Propiedades del Esquema
3.2.1 Motivación y contexto metodológico
La resolución analítica de ecuaciones diferenciales estocásticas con retraso (SDDE) es generalmente imposible debido a la naturaleza no markoviana del proceso y la dependencia funcional del pasado (Kloeden & Platen, 1992). Por lo tanto, es necesario recurrir a métodos numéricos de discretización para obtener soluciones aproximadas que permitan tanto el análisis teórico como la implementación computacional (Mao, 2007).
El esquema de Euler-Maruyama representa el método más fundamental y ampliamente utilizado para la discretización de SDEs y SDDEs, debido a su simplicidad conceptual, facilidad de implementación y propiedades de convergencia bien establecidas (Platen, 1999).
3.2.2 Esquema de Euler-Maruyama para SDDE
Consideremos la SDDE definida en la sección 3.1:
\[ dX(t) = \mu(t, X(t), X(t-\tau)) \, dt + \sigma(t, X(t)) \, dW(t), \quad t \geq 0 \tag{3.1} \]
donde \(\mu(t, x, y) = \alpha + \beta x + \gamma y + \boldsymbol{\delta}^T \mathbf{C}(t)\) y \(\sigma(t, x) = \sigma x\).
Para discretizar esta ecuación, particionamos el intervalo temporal \([0, T]\) en \(N\) subintervalos de igual longitud \(\Delta t = T/N\), definiendo los puntos de malla \(t_n = n \Delta t\) para \(n = 0, 1, \dots, N\). Denotamos \(X_n \approx X(t_n)\) como la aproximación numérica del proceso en el tiempo \(t_n\).
El esquema de Euler-Maruyama para la SDDE (3.1) se define recursivamente como:
\[ X_{n+1} = X_n + \mu(t_n, X_n, X_{n-k}) \Delta t + \sigma(t_n, X_n) \Delta W_n \]
donde: - \(k = \tau / \Delta t\) es el número de pasos de retraso (asumiendo que \(\tau\) es múltiplo entero de \(\Delta t\)), - \(\Delta W_n = W(t_{n+1}) - W(t_n) \sim \mathcal{N}(0, \Delta t)\) es el incremento del proceso de Wiener, - \(X_{n-k} \approx X(t_n - \tau)\) es la aproximación del estado retardado.
Para nuestro modelo específico con \(\mu(t, x, y) = \alpha + \beta x + \gamma y + \boldsymbol{\delta}^T \mathbf{C}(t)\) y \(\sigma(t, x) = \sigma x\), el esquema se particulariza como:
\[ X_{n+1} = X_n + \left[ \alpha + \beta X_n + \gamma X_{n-k} + \boldsymbol{\delta}^T \mathbf{C}_n \right] \Delta t + \sigma X_n \Delta W_n \tag{3.10} \]
donde \(\mathbf{C}_n = \mathbf{C}(t_n)\) representa el vector de forzantes climáticas en el tiempo discreto \(t_n\).
3.2.3 Representación alternativa en diferencias
En aplicaciones prácticas con datos observados, es conveniente reorganizar la ecuación (3.10) en términos del diferencial discreto \(\Delta X_n = X_{n+1} - X_n\):
\[ \Delta X_n = \alpha \Delta t + \beta X_n \Delta t + \gamma X_{n-k} \Delta t + \boldsymbol{\delta}^T \mathbf{C}_n \Delta t + \sigma X_n \Delta W_n \tag{3.11} \]
Esta formulación es particularmente útil para la estimación de parámetros mediante métodos de regresión, ya que permite expresar el cambio en el estado como una función lineal de las variables explicativas (Shoji & Ozaki, 1998).
3.2.4 Análisis de convergencia
El análisis de convergencia de esquemas numéricos para SDEs y SDDEs distingue entre dos tipos fundamentales de convergencia: convergencia fuerte y convergencia débil (Kloeden & Platen, 1992).
3.2.4.1 Convergencia fuerte
Un esquema numérico tiene convergencia fuerte de orden \(\gamma\) si existe una constante \(C > 0\) independiente de \(\Delta t\) tal que:
\[ \mathbb{E}\left[ |X(T) - X_N| \right] \leq C (\Delta t)^\gamma \]
donde \(X(T)\) es la solución exacta en tiempo \(T\) y \(X_N\) es la aproximación numérica.
Teorema 3.1 (Convergencia fuerte de Euler-Maruyama). Bajo las condiciones de Lipschitz global y crecimiento lineal para \(\mu\) y \(\sigma\), el esquema de Euler-Maruyama (3.10) tiene convergencia fuerte de orden \(\gamma = 1/2\) (Kloeden & Platen, 1992).
Demostración: Ver (Kloeden & Platen, 1992), Teorema 10.2.2, o Anexo D para una demostración adaptada a SDDEs.
Para nuestro modelo lineal específico, la condición de Lipschitz se satisface trivialmente, garantizando la convergencia fuerte de orden \(1/2\).
3.2.4.2 Convergencia débil
Un esquema numérico tiene convergencia débil de orden \(\beta\) si para toda función \(g \in C_P^{2(\beta+1)}(\mathbb{R}, \mathbb{R})\) (funciones \(2(\beta+1)\) veces continuamente diferenciables con derivadas polinomialmente acotadas), existe una constante \(C_g > 0\) tal que:
\[ \left| \mathbb{E}[g(X(T))] - \mathbb{E}[g(X_N)] \right| \leq C_g (\Delta t)^\beta \tag{3.12} \]
Teorema 3.2 (Convergencia débil de Euler-Maruyama). Bajo las mismas condiciones del Teorema 3.1, el esquema de Euler-Maruyama tiene convergencia débil de orden \(\beta = 1\) (Milstein, 1995).
Este resultado es particularmente relevante para aplicaciones en finanzas y ciencias naturales, donde el interés radica en calcular esperanzas de funcionales del proceso (precios de opciones, probabilidades de eventos, etc.) (Glasserman, 2003).
3.2.5 Estabilidad numérica
La estabilidad numérica de un esquema de discretización se refiere a su capacidad para mantener acotadas las soluciones numéricas cuando las soluciones exactas son acotadas (Higham et al., 2002). Para el esquema de Euler-Maruyama aplicado a SDEs lineales, se pueden derivar condiciones explícitas de estabilidad.
Consideremos el caso simplificado sin retraso (\(\gamma = 0\)) y sin forzantes climáticas (\(\boldsymbol{\delta} = \mathbf{0}\)):
\[ X_{n+1} = X_n + (\alpha + \beta X_n) \Delta t + \sigma X_n \Delta W_n \tag{3.13} \]
Reorganizando términos:
\[ X_{n+1} = \left[ 1 + \beta \Delta t + \sigma \Delta W_n \right] X_n + \alpha \Delta t \tag{3.14} \]
Definición 3.1 (Estabilidad MS). Un esquema numérico es estable en media cuadrática (MS) si existe una constante \(C > 0\) tal que:
\[ \sup_{n \geq 0} \mathbb{E}[|X_n|^2] \leq C \mathbb{E}[|X_0|^2] \tag{3.15} \]
Teorema 3.3 (Condición de estabilidad MS para Euler-Maruyama). El esquema (3.14) es estable en media cuadrática si y solo si:
\[ 2\beta + \sigma^2 < 0 \tag{3.16} \]
Demostración: Ver (Higham et al., 2002), Teorema 2.1.
Esta condición establece que la estabilidad del esquema depende del balance entre el coeficiente de retroalimentación \(\beta\) y la intensidad del ruido \(\sigma\). Intuitivamente, el ruido puede tener un efecto estabilizador o desestabilizador dependiendo de la magnitud relativa de estos parámetros (Mao, 2007).
3.2.6 Error de truncamiento local
El error de truncamiento local en el paso \(n\) se define como la diferencia entre la solución exacta y la aproximación numérica, asumiendo que el estado anterior es exacto:
\[ R_n = X(t_{n+1}) - X(t_n) - \mu(t_n, X(t_n), X(t_n-\tau)) \Delta t - \sigma(t_n, X(t_n)) \Delta W_n \tag{3.17} \]
Para el esquema de Euler-Maruyama aplicado a SDEs con coeficientes suficientemente suaves, se puede demostrar que:
\[ \mathbb{E}[|R_n|^2]^{1/2} = \mathcal{O}((\Delta t)^{3/2}) \tag{3.18} \]
Este resultado justifica el orden de convergencia fuerte \(\gamma = 1/2\), ya que el error acumulado sobre \(N = T/\Delta t\) pasos es de orden \(\mathcal{O}((\Delta t)^{1/2})\) (Kloeden & Platen, 1992).
3.2.7 Implementación computacional
La implementación del esquema de Euler-Maruyama para SDDEs requiere especial atención al manejo del historial pasado. El Algoritmo 3.1 presenta una implementación eficiente:
Algoritmo 3.1: Esquema de Euler-Maruyama para SDDE
Entrada: \(\tau, \Delta t, T, X_0, \{C_n\}_{n=0}^{N-1},\alpha, \beta, \gamma, \delta, \sigma\).
Salida: \(\{X_n\}_{n=0}^{N}\) aproximación numérica.
\(N\longleftarrow \dfrac{T}{\Delta t}\) \(k \longleftarrow \dfrac{\tau}{\Delta t}\) Inicializar vector \(X[0..N]\) \(X[0..k]\longleftarrow \phi(t_0..t_k)\) // Historial inicial Para \(n = k\) hasta \(N-1\) hacer: a. \(\Delta W \longleftarrow \sqrt{\Delta t}\cdot \mathcal{N}(0, 1)\) b. \(\mu \longleftarrow \alpha + \beta\cdot X[n] + \gamma\cdot X[n-k] + \mathbf{\delta}^T C[n]\) c. \(X[n+1] \longleftarrow X[n] + \mu\cdot \Delta t + \sigma\cdot X[n] \Delta W\) Retornar \(\{X_n\}_{n=0}^{N}\)
La complejidad computacional del algoritmo es \(\mathcal{O}(N)\), lo que lo hace eficiente incluso para simulaciones con horizontes temporales largos (andersson2017strong?).
3.2.8 Consideraciones prácticas para \(\Delta t = 1\) día
En nuestra aplicación específica, el paso temporal de discretización es \(\Delta t = 1\) día, correspondiente a la frecuencia de muestreo de los datos PhenoCam. Esta elección tiene implicaciones importantes:
Magnitud del error numérico: Con \(\Delta t = 1\), el error de truncamiento local puede ser significativo comparado con aplicaciones donde \(\Delta t \ll 1\). Sin embargo, dado que los datos observados también tienen una resolución diaria, el error numérico se vuelve comparable al error de observación (Shoji & Ozaki, 1998).
Aproximación del incremento de Wiener: Para \(\Delta t = 1\), tenemos \(\Delta W_n \sim \mathcal{N}(0, 1)\), lo que simplifica la generación de números aleatorios en la simulación.
Interpretación de parámetros: Los parámetros estimados (\(\alpha\), \(\beta\), \(\gamma\), \(\boldsymbol{\delta}\), \(\sigma\)) están en escala diaria y deben interpretarse como tasas de cambio por día, no como tasas instantáneas (Campbell et al., 1997).
Consistencia dimensional: Con \(\Delta t = 1\), las ecuaciones (3.9) y (3.10) se simplifican al omitir explícitamente el factor \(\Delta t\), lo que facilita la interpretación y la estimación de parámetros.
3.3 Formulación Estadística como Modelo de Regresión
3.3.1 Motivación y contexto estadístico
La discretización del modelo SDDE mediante el esquema de Euler-Maruyama, presentada en la sección 3.2, proporciona una aproximación numérica que facilita la transición desde un marco teórico continuo a un problema de inferencia estadística discreta (Shoji & Ozaki, 1998). Esta formulación permite estimar los parámetros del modelo utilizando técnicas estándar de regresión lineal, lo cual es computacionalmente eficiente y estadísticamente bien fundamentado (Campbell et al., 1997).
El enfoque de regresión lineal para modelos estocásticos discretizados ha sido ampliamente utilizado en finanzas, ecología y ciencias ambientales, donde la combinación de dinámica determinista y ruido estocástico es inherente a los sistemas naturales (Kloeden & Platen, 1992).
3.3.2 Modelo de regresión lineal derivado del SDDE
Partiendo del esquema de Euler-Maruyama para nuestro modelo SDDE específico (ecuación 3.9), y considerando un paso temporal \(\Delta t = 1\) día, podemos reorganizar la ecuación para expresar el cambio en el estado como una función lineal de las variables explicativas:
\[ \Delta X_n = \alpha + \beta X_n + \gamma X_{n-k} + \delta_1 \text{Rad}_n + \delta_2 \text{Temp}_n + \delta_3 \text{Prec}_n + \epsilon_n \tag{3.19} \]
donde: - \(\Delta X_n = X_{n+1} - X_n\) es el cambio observado en el índice GCC entre el día \(n\) y \(n+1\), - \(X_n\) es el valor del GCC en el día \(n\), - \(X_{n-k}\) es el valor del GCC en el día \(n-k\), donde \(k = \tau\) representa el número de días de retraso, - \(\text{Rad}_n\), \(\text{Temp}_n\), \(\text{Prec}_n\) son las variables climáticas observadas en el día \(n\), - \(\epsilon_n = \sigma X_n \Delta W_n\) es el término de error estocástico.
Esta formulación constituye un modelo de regresión lineal múltiple donde la variable dependiente es \(\Delta X_n\) y las variables independientes son \(X_n\), \(X_{n-k}\), y las tres variables climáticas.
3.3.3 Notación matricial y vector de parámetros
En notación matricial compacta, el modelo (3.19) se expresa como:
\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\epsilon} \tag{3.20} \]
donde: - \(\mathbf{Y} \in \mathbb{R}^N\) es el vector de observaciones de la variable dependiente (\(\Delta X_n\)), - \(\mathbf{X} \in \mathbb{R}^{N \times p}\) es la matriz de diseño con \(p = 6\) columnas, - \(\boldsymbol{\theta} \in \mathbb{R}^p\) es el vector de parámetros desconocidos, - \(\boldsymbol{\epsilon} \in \mathbb{R}^N\) es el vector de errores aleatorios.
La matriz de diseño \(\mathbf{X}\) tiene la siguiente estructura:
\[ \mathbf{X} = \begin{bmatrix} 1 & X_1 & X_{1-k} & \text{Rad}_1 & \text{Temp}_1 & \text{Prec}_1 \\ 1 & X_2 & X_{2-k} & \text{Rad}_2 & \text{Temp}_2 & \text{Prec}_2 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & X_N & X_{N-k} & \text{Rad}_N & \text{Temp}_N & \text{Prec}_N \end{bmatrix} \tag{3.21} \]
y el vector de parámetros es:
\[ \boldsymbol{\theta} = \begin{bmatrix} \alpha & \beta & \gamma & \delta_1 & \delta_2 & \delta_3 \end{bmatrix}^T \tag{3.22} \]
3.3.4 Relación estructural con modelos AR(p)
El modelo de regresión (3.19) puede reinterpretarse como un modelo autorregresivo con estructura esparsa. Reorganizando la ecuación (3.19):
\[ X_{n+1} = \alpha + (1 + \beta) X_n + \gamma X_{n-k} + \delta_1 \text{Rad}_n + \delta_2 \text{Temp}_n + \delta_3 \text{Prec}_n + \epsilon_n \tag{3.23} \]
Comparando con un modelo AR(p) general:
\[ X_{n+1} = c + \sum_{i=1}^p \phi_i X_{n+1-i} + \varepsilon_n \tag{3.24} \]
podemos establecer la siguiente correspondencia estructural:
\[ \phi_i = \begin{cases} 1 + \beta & \text{si } i = 1 \\ \gamma & \text{si } i = k \\ 0 & \text{en otro caso} \end{cases} \tag{3.25} \]
Esta relación demuestra que nuestro modelo SDDE discretizado es equivalente a un modelo AR(p) con restricciones estructurales, donde solo dos coeficientes autorregresivos son no nulos: el correspondiente al rezago inmediato (\(i=1\)) y el correspondiente al retraso biológico (\(i=k=\tau\)). Esta estructura esparsa reduce drásticamente la dimensionalidad del problema de \(p+1\) parámetros a solo 6 parámetros, imponiendo una hipótesis biológicamente significativa sobre la memoria temporal del sistema (Vapnik, 1998).
3.3.5 Supuestos estadísticos del modelo
Para que los estimadores de mínimos cuadrados ordinarios (OLS) tengan propiedades óptimas, se requieren los siguientes supuestos clásicos del modelo lineal (greene2018econometric?):
Supuesto 1: Linealidad en parámetros El modelo es lineal en los parámetros \(\boldsymbol{\theta}\), lo cual se satisface por construcción.
Supuesto 2: Rango completo de la matriz de diseño La matriz \(\mathbf{X}\) tiene rango completo (\(\text{rank}(\mathbf{X}) = p\)), lo que garantiza que \(\mathbf{X}^T\mathbf{X}\) sea invertible. Este supuesto requiere que no exista multicolinealidad perfecta entre las variables explicativas.
Supuesto 3: Esperanza condicional nula \[ \mathbb{E}[\epsilon_n | \mathbf{X}] = 0 \tag{3.26} \] Este supuesto implica que los errores no están correlacionados con las variables explicativas, lo cual es crucial para la insesgadez de los estimadores.
Supuesto 4: Homocedasticidad y no autocorrelación \[ \text{Var}(\boldsymbol{\epsilon} | \mathbf{X}) = \sigma^2 \mathbf{I}_N \tag{3.27} \] donde \(\mathbf{I}_N\) es la matriz identidad de dimensión \(N\). Este supuesto implica que los errores tienen varianza constante y no están autocorrelacionados.
Supuesto 5: Normalidad de los errores (para inferencia exacta) \[ \boldsymbol{\epsilon} | \mathbf{X} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_N) \tag{3.28} \] Este supuesto es necesario para realizar pruebas de hipótesis y construir intervalos de confianza exactos, aunque no es requerido para la consistencia de los estimadores.
3.3.6 Estimación de parámetros por mínimos cuadrados ordinarios
Bajo los supuestos anteriores, el estimador de mínimos cuadrados ordinarios (OLS) de \(\boldsymbol{\theta}\) se define como:
\[ \hat{\boldsymbol{\theta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} \tag{3.29} \]
Este estimador tiene las siguientes propiedades fundamentales (greene2018econometric?):
- Insesgadez: \(\mathbb{E}[\hat{\boldsymbol{\theta}}] = \boldsymbol{\theta}\)
- Consistencia: \(\hat{\boldsymbol{\theta}} \xrightarrow{p} \boldsymbol{\theta}\) cuando \(N \to \infty\)
- Eficiencia: Entre todos los estimadores lineales insesgados, \(\hat{\boldsymbol{\theta}}\) tiene la menor varianza (teorema de Gauss-Markov)
La matriz de varianza-covarianza del estimador es:
\[ \text{Var}(\hat{\boldsymbol{\theta}}) = \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1} \tag{3.30} \]
donde \(\sigma^2\) se estima mediante:
\[ \hat{\sigma}^2 = \frac{1}{N - p} \sum_{n=1}^N \hat{\epsilon}_n^2 = \frac{\|\mathbf{Y} - \mathbf{X} \hat{\boldsymbol{\theta}}\|^2}{N - p} \tag{3.31} \]
3.3.7 Identificabilidad del modelo
La identificabilidad de un modelo estadístico se refiere a la capacidad de estimar de manera única sus parámetros a partir de los datos observados (rothenberg1971identification?). Para nuestro modelo, la identificabilidad depende de dos factores principales:
Identificabilidad estructural El modelo es estructuralmente identificable si la parametrización es única, es decir, si diferentes valores de \(\boldsymbol{\theta}\) producen distribuciones de probabilidad diferentes para los datos observados. En nuestro caso, la linealidad del modelo y la independencia funcional de las variables explicativas garantizan la identificabilidad estructural.
Identificabilidad práctica La identificabilidad práctica depende de la calidad y cantidad de los datos observados. Específicamente, se requiere que:
La matriz \(\mathbf{X}^T \mathbf{X}\) sea no singular (rango completo)
Las variables climáticas no estén perfectamente colineales con los términos de retraso
Exista suficiente variabilidad en los datos para distinguir los efectos individuales de cada parámetro
En la práctica, la identificabilidad se evalúa mediante el factor de inflación de varianza (VIF) y el número de condición de la matriz \(\mathbf{X}^T \mathbf{X}\). Un número de condición elevado (> 1000) indica problemas de identificabilidad práctica debido a multicolinealidad (belsley1980regression?).
3.3.8 Incorporación de variables exógenas climáticas
La inclusión de variables climáticas exógenas en el modelo responde a dos objetivos fundamentales:
Control de factores de confusión: Las variables climáticas pueden estar correlacionadas tanto con el estado actual del sistema como con su estado pasado, por lo que su inclusión evita sesgos en la estimación de los coeficientes de retroalimentación \(\beta\) y \(\gamma\).
Mejora de la capacidad predictiva: Las variables climáticas contienen información adicional sobre los forzantes ambientales que afectan directamente el crecimiento vegetal, lo que mejora la precisión de las predicciones del modelo.
La especificación funcional lineal adoptada para las variables climáticas asume que sus efectos son aditivos y constantes en el tiempo. Esta suposición puede ser relajada en extensiones futuras del modelo mediante la inclusión de términos no lineales o interacciones (hastie2009elements?).
3.3.9 Diagnóstico de supuestos y validación del modelo
La validez de las inferencias estadísticas basadas en el modelo de regresión depende de la verificación de los supuestos establecidos en la sección 3.3.5. Los principales métodos de diagnóstico incluyen:
- Gráfico de residuos vs valores ajustados: Para detectar heterocedasticidad y no linealidad
- Test de Durbin-Watson: Para detectar autocorrelación serial en los residuos
- Test de Breusch-Pagan: Para detectar heterocedasticidad
- QQ-plot de residuos: Para evaluar la normalidad de los errores
- Análisis de influencia: Para identificar observaciones atípicas que puedan afectar desproporcionadamente los resultados
Estos diagnósticos son esenciales para garantizar la robustez de las conclusiones derivadas del modelo y para identificar posibles violaciones de los supuestos que requieran modificaciones en la especificación del modelo (fox2015applied?).
3.4 Identificabilidad y Estimación de Parámetros
3.4.1 Marco teórico de identificabilidad
La identificabilidad constituye una propiedad fundamental en la inferencia estadística, ya que garantiza que los parámetros de un modelo puedan ser estimados de manera única a partir de los datos observados (rothenberg1971identification?). En el contexto de modelos estocásticos discretizados como el nuestro, la identificabilidad se descompone en dos componentes complementarios: identificabilidad estructural e identificabilidad práctica.
La identificabilidad estructural se refiere a la unicidad de la parametrización del modelo desde una perspectiva puramente matemática, independientemente de los datos disponibles. Por otro lado, la identificabilidad práctica depende de la calidad, cantidad y estructura de los datos observados, y determina si es posible estimar los parámetros con precisión finita (bellman1970structural?).
3.4.2 Identificabilidad estructural del modelo SDDE
Consideremos el modelo de regresión lineal derivado del SDDE (ecuación 3.19):
\[ \Delta X_n = \alpha + \beta X_n + \gamma X_{n-k} + \delta_1 \text{Rad}_n + \delta_2 \text{Temp}_n + \delta_3 \text{Prec}_n + \epsilon_n \tag{3.32} \]
El modelo es estructuralmente identificable si la aplicación \(\boldsymbol{\theta} \mapsto \mathbb{P}_{\boldsymbol{\theta}}\) es inyectiva, donde \(\mathbb{P}_{\boldsymbol{\theta}}\) denota la distribución de probabilidad de los datos bajo el parámetro \(\boldsymbol{\theta}\) (walter1982global?).
Teorema 3.4 (Identificabilidad estructural). El modelo (3.32) es estructuralmente identificable si y solo si las funciones \(X_n\), \(X_{n-k}\), \(\text{Rad}_n\), \(\text{Temp}_n\), y \(\text{Prec}_n\) son linealmente independientes en el espacio de funciones generadas por el proceso estocástico.
Demostración. La linealidad del modelo implica que si existen dos vectores de parámetros \(\boldsymbol{\theta}_1 \neq \boldsymbol{\theta}_2\) tales que \(\mathbb{P}_{\boldsymbol{\theta}_1} = \mathbb{P}_{\boldsymbol{\theta}_2}\), entonces:
\[ (\boldsymbol{\theta}_1 - \boldsymbol{\theta}_2)^T \mathbf{x}_n = 0 \quad \text{casi seguramente para todo } n \tag{3.33} \]
donde \(\mathbf{x}_n = [1, X_n, X_{n-k}, \text{Rad}_n, \text{Temp}_n, \text{Prec}_n]^T\). Esto ocurre si y solo si las componentes de \(\mathbf{x}_n\) son linealmente dependientes, lo cual contradice la hipótesis del teorema. ∎
En la práctica, la independencia lineal se satisface siempre que: 1. Las variables climáticas no son funciones deterministas del estado del sistema (\(X_n\) o \(X_{n-k}\)) 2. El retraso \(\tau\) es tal que \(X_n\) y \(X_{n-\tau}\) no son perfectamente correlacionados 3. Existe suficiente variabilidad en los datos para distinguir los efectos individuales
3.4.3 Identificabilidad práctica y condiciones de rango
La identificabilidad práctica se evalúa mediante el rango de la matriz de información de Fisher o, equivalentemente, el rango de la matriz de diseño \(\mathbf{X}\) (greene2018econometric?).
Definición 3.2 (Condición de rango completo). El modelo es prácticamente identificable si la matriz \(\mathbf{X}^T \mathbf{X}\) es no singular, es decir, si \(\text{rank}(\mathbf{X}) = p = 6\).
Esta condición es equivalente a requerir que el determinante de Gram sea no nulo:
\[ \det(\mathbf{X}^T \mathbf{X}) > 0 \tag{3.34} \]
En términos geométricos, esto significa que los vectores columna de \(\mathbf{X}\) son linealmente independientes en \(\mathbb{R}^N\).
Proposición 3.1 (Condiciones suficientes para identificabilidad práctica). El modelo (3.32) es prácticamente identificable si:
- \(N > p\) (más observaciones que parámetros)
- No existe multicolinealidad perfecta entre las variables explicativas
- Las variables climáticas presentan variabilidad suficiente y no están perfectamente correlacionadas con los términos de estado
La violación de estas condiciones conduce a problemas de multicolinealidad, que se manifiestan en varianzas infladas de los estimadores y dificultades numéricas en la inversión de \(\mathbf{X}^T \mathbf{X}\) (belsley1980regression?).
3.4.4 Estimación por mínimos cuadrados ordinarios (OLS)
Bajo las condiciones de identificabilidad establecidas, el estimador de mínimos cuadrados ordinarios (OLS) se define como la solución al problema de optimización:
\[ \hat{\boldsymbol{\theta}} = \arg\min_{\boldsymbol{\theta} \in \mathbb{R}^p} \|\mathbf{Y} - \mathbf{X} \boldsymbol{\theta}\|^2 \tag{3.35} \]
La solución analítica a este problema viene dada por la ecuación normal:
\[ \hat{\boldsymbol{\theta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} \tag{3.36} \]
Este estimador posee propiedades óptimas bajo los supuestos del modelo lineal clásico:
- Insesgadez: \(\mathbb{E}[\hat{\boldsymbol{\theta}}] = \boldsymbol{\theta}\)
- Consistencia: \(\hat{\boldsymbol{\theta}} \xrightarrow{p} \boldsymbol{\theta}\) cuando \(N \to \infty\)
- Eficiencia: \(\text{Var}(\hat{\boldsymbol{\theta}}) \leq \text{Var}(\tilde{\boldsymbol{\theta}})\) para cualquier otro estimador lineal insesgado \(\tilde{\boldsymbol{\theta}}\) (teorema de Gauss-Markov)
La matriz de varianza-covarianza del estimador es:
\[ \text{Var}(\hat{\boldsymbol{\theta}}) = \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1} \tag{3.37} \]
donde \(\sigma^2 = \text{Var}(\epsilon_n)\) es la varianza del error.
3.4.5 Estimación de la varianza del error
La varianza del error \(\sigma^2\) se estima mediante el estimador insesgado:
\[ \hat{\sigma}^2 = \frac{1}{N - p} \sum_{n=1}^N \hat{\epsilon}_n^2 = \frac{\|\mathbf{Y} - \mathbf{X} \hat{\boldsymbol{\theta}}\|^2}{N - p} \tag{3.38} \]
donde \(\hat{\epsilon}_n = Y_n - \mathbf{x}_n^T \hat{\boldsymbol{\theta}}\) son los residuos estimados.
Este estimador es insesgado bajo el supuesto de homocedasticidad y tiene la propiedad de que \((N - p) \hat{\sigma}^2 / \sigma^2 \sim \chi^2_{N - p}\) bajo normalidad de los errores (casella2002statistical?).
3.4.6 Relación con la intensidad del ruido estocástico
En el contexto del modelo SDDE original, la varianza del error \(\sigma^2\) está relacionada con la intensidad del ruido estocástico \(\sigma_{\text{SDDE}}\) del proceso continuo. Recordando que el término de error en la discretización es \(\epsilon_n = \sigma_{\text{SDDE}} X_n \Delta W_n\), y que \(\Delta W_n \sim \mathcal{N}(0, \Delta t)\), tenemos:
\[ \text{Var}(\epsilon_n | X_n) = \sigma_{\text{SDDE}}^2 X_n^2 \Delta t \tag{3.39} \]
Sin embargo, en nuestra formulación de regresión asumimos homocedasticidad (\(\text{Var}(\epsilon_n) = \sigma^2\)), lo cual constituye una aproximación válida cuando la variabilidad relativa de \(X_n\) es pequeña o cuando se trabaja con tasas de crecimiento relativas en lugar de absolutas (Shoji & Ozaki, 1998).
Para recuperar una estimación de la intensidad del ruido del proceso continuo, podemos utilizar:
\[ \hat{\sigma}_{\text{SDDE}} = \sqrt{\frac{\hat{\sigma}^2}{\frac{1}{N} \sum_{n=1}^N X_n^2 \Delta t}} \tag{3.40} \]
Esta estimación es consistente bajo condiciones regulares y permite interpretar el parámetro \(\sigma_{\text{SDDE}}\) en términos del modelo continuo original.
3.4.7 Diagnóstico de identificabilidad práctica
La identificabilidad práctica se diagnostica mediante métricas numéricas que cuantifican el grado de multicolinealidad en la matriz de diseño:
3.4.7.1 Factor de inflación de varianza (VIF)
Para cada variable explicativa \(j\), el VIF se define como:
\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \tag{3.41} \]
donde \(R_j^2\) es el coeficiente de determinación obtenido al regresar la variable \(j\) contra todas las demás variables explicativas. Valores de \(\text{VIF}_j > 10\) indican problemas severos de multicolinealidad (marquardt1970generalized?).
3.4.7.2 Número de condición
El número de condición de la matriz \(\mathbf{X}^T \mathbf{X}\) se define como:
\[ \kappa(\mathbf{X}^T \mathbf{X}) = \sqrt{\frac{\lambda_{\max}}{\lambda_{\min}}} \tag{3.42} \]
donde \(\lambda_{\max}\) y \(\lambda_{\min}\) son los autovalores máximo y mínimo, respectivamente. Según (belsley1980regression?):
- \(\kappa < 10\): Poca multicolinealidad
- \(10 \leq \kappa < 30\): Multicolinealidad moderada
- \(\kappa \geq 30\): Multicolinealidad severa
Estos diagnósticos son esenciales para evaluar la fiabilidad de las estimaciones de parámetros y para identificar posibles reformulaciones del modelo que mejoren su identificabilidad práctica.