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:

\[ \boxed{ 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{4.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{4.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{4.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 (4.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{4.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{4.5} \]

Teorema 4.1 (Existencia y unicidad). Bajo las condiciones (4.4) y (4.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 (4.1) y la condición inicial (4.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:

  1. 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|)\).

  2. 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 (4.1) puede descomponerse en cuatro componentes:

\[ \mu(t, X(t), X(t-\tau)) = \underbrace{\alpha}_{\text{Crecimiento basal}} + \underbrace{\beta X(t)}_{\text{Retroalimentación actual}} + \underbrace{\gamma X(t-\tau)}_{\text{Memoria fisiológica}} + \underbrace{\boldsymbol{\delta}^T \mathbf{C}(t)}_{\text{Forzantes climáticas}} \tag{4.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 (4.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{4.7} \]

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 (4.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 (4.7) 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{4.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{4.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.7} \]

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.7) se define recursivamente como:

\[ \boxed{ X_{n+1} = X_n + \mu(t_n, X_n, X_{n-k}) \Delta t + \sigma(t_n, X_n) \Delta W_n } \tag{3.8} \]

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:

\[ \boxed{ 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.9} \]

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.9) en términos del diferencial discreto \(\Delta X_n = X_{n+1} - X_n\):

\[ \boxed{ \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.10} \]

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 \tag{3.11} \]

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.8) 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:

\[ \boxed{ 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.

  1. \(N\longleftarrow \dfrac{T}{\Delta t}\)
  2. \(k \longleftarrow \dfrac{\tau}{\Delta t}\)
  3. Inicializar vector \(X[0..N]\)
  4. \(X[0..k]\longleftarrow \phi(t_0..t_k)\) // Historial inicial
  5. Para \(n = k\) hasta \(N-1\) hacer:
    1. \(\Delta W \longleftarrow \sqrt{\Delta t}\cdot \mathcal{N}(0, 1)\)
    2. \(\mu \longleftarrow \alpha + \beta\cdot X[n] + \gamma\cdot X[n-k] + \mathbf{\delta}^T C[n]\)
    3. \(X[n+1] \longleftarrow X[n] + \mu\cdot \Delta t + \sigma\cdot X[n] \Delta W\)
  6. 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:

  1. 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).

  2. 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.

  3. 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).

  4. 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.

Campbell, J. Y., Lo, A. W., & MacKinlay, A. C. (1997). The Econometrics of Financial Markets. Princeton University Press.
Diepenbrock, W. (2000). Yield formation in field-grown Brassica oleracea L. var. botrytis: A review. Journal of Horticultural Science & Biotechnology, 75(4), 395-408. https://doi.org/10.1080/14620316.2000.11511267
Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering (Vol. 53). Springer. https://doi.org/10.1007/978-0-387-21617-1
Higham, D. J., Mao, X., & Stuart, A. M. (2002). Strong convergence of Euler-type methods for nonlinear stochastic differential equations. SIAM Journal on Numerical Analysis, 40(3), 1041-1063. https://doi.org/10.1137/S0036142901389530
Itô, K. (1951). On stochastic differential equations. Memoirs of the American Mathematical Society, 4, 1-51. https://doi.org/10.1090/memo/0004
Kloeden, P. E., & Platen, E. (1992). Numerical Solution of Stochastic Differential Equations (Vol. 23). Springer-Verlag. https://doi.org/10.1007/978-3-662-12616-5
Mao, X. (2007). Stochastic Differential Equations and Applications (2.ª ed.). Horwood Publishing.
Milstein, G. N. (1995). Numerical Integration of Stochastic Differential Equations (Vol. 313). Kluwer Academic Publishers. https://doi.org/10.1007/978-94-015-8455-5
Mohammed, S.-E. A. (1984). Stochastic Functional Differential Equations. 99.
Monteith, J. L., & Moss, C. J. (1977). Climate and the efficiency of crop production in Britain. Philosophical Transactions of the Royal Society B: Biological Sciences, 281(980), 277-294. https://doi.org/10.1098/rstb.1977.0140
Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications (6.ª ed.). Springer. https://doi.org/10.1007/978-3-642-14394-6
Platen, E. (1999). An introduction to numerical methods for stochastic differential equations. Acta Numerica, 8, 197-246. https://doi.org/10.1017/S0962492900002941
Revuz, D., & Yor, M. (1999). Continuous Martingales and Brownian Motion (3.ª ed., Vol. 293). Springer.
Shoji, I., & Ozaki, T. (1998). Estimation for nonlinear stochastic differential equations by a local linearization method. Stochastic Analysis and Applications, 16(4), 733-752. https://doi.org/10.1080/07362999808809559
Taiz, L., & Zeiger, E. (2010). Plant Physiology (5.ª ed.). Sinauer Associates.
Thomas, H., & Ougham, H. (2014). The stay-green trait. Journal of Experimental Botany, 65(14), 3889-3900. https://doi.org/10.1093/jxb/eru033
Volterra, V. (1931). Théorie Mathématique de la Lutte pour la Vie. Gauthier-Villars.