Resulta que podemos obtener una solución analítica exacta.
Porque $(x_t, x_{t+1}, x_{t+2})$ tiene una distribución Normal (trivariada) (como se demostrará más adelante), cualquier combinación lineal también es Normal, como por ejemplo $(X,Y)=(x_{t+1}-x_t, x_{t+2}-x_t).$ Es evidente que tiene una expectativa nula, por lo que su distribución vendrá determinada por su matriz de varianzas-covarianzas,
$$\operatorname{Var}((X,Y)) = \Sigma_{X,Y} = \pmatrix{\operatorname{Var}(x_{t+1}-x_t) & \operatorname{Cov}(x_{t+1}-x_t,\ x_{t+2}-x_t\\\operatorname{Cov}(x_{t+2}-x_t,\ x_{t+1}-x_t) & \operatorname{Var}(x_{t+2}-x_t)}.$$
Utilice la definición recursiva de la serie para calcular estos valores, empezando por
$$\begin{aligned} x_{t+1} - x_t &= (b-1)x_t + u_{t+1};\\ x_{t+2}-x_t &= b\,x_{t+1}+u_{t+2} - x_t = b(b\,x_t + u_{t+1}) + u_{t+2} - x_t \\&= (b^2-1)x_t + b\,u_{t+1}+u_{t+2}. \end{aligned}$$
Todo es en términos de combinaciones lineales de los independiente Variables normales $(x_t, u_{t+1}, u_{t+2}).$ Esto hace que $(X,Y)$ Normal multivariante, como se ha afirmado anteriormente. Pasemos a calcular $\Sigma_{X,Y}.$
Presumiblemente la serie $(x_t)$ es estacionario, lo que implica $\operatorname{Var}(x_t) = \sigma^2$ es el mismo para todos $t.$ Así
$$\sigma^2 = \operatorname{Var}(x_{t}) = \operatorname{Var}(bx_{t-1} + u_{t}) = b^2\sigma^2 + 1,$$
con la solución
$$\operatorname{Var}(x_{t}) = \sigma^2 = \frac{1}{1-b^2}.$$
El resto del cálculo de la varianza encaja a la perfección:
$$\operatorname{Var}(x_{t+1}-x_t) = \operatorname{Var}((b-1)x_t + u_{t+1}) = (b-1)^2\sigma^2 + 1 = \frac{2}{1+b}; $$ $$\operatorname{Var}(x_{t+2}-x_t) = \operatorname{Var}((b^2-1)x_t + bu_{t+1} + u_{t+2}) = (b^2-1)^2\sigma^2 + b^2 + 1 = 2$$ $$\operatorname{Cov}(x_{t+1}-x_t, x_{t+2}-x_t) = (b-1)(b^2-1)/\sigma^2 + b = 1.$$
Eso es, la matriz de varianza-covarianza de $(X,Y)$ es simplemente
$$\Sigma_{X,Y} = \pmatrix{\frac{2}{1+b} & 1 \\ 1 & 2}.$$
Por lo tanto, su coeficiente de correlación es
$$\rho = \frac{1}{\sqrt{2/(1+b)\times 2}} = \frac{\sqrt{1+b}}{2}.$$
El método de solución en https://stats.stackexchange.com/a/139090/919 para hallar la densidad del máximo de dos variables normales correlacionadas con igual varianzas, puede extenderse directamente a una distribución Normal bivariante arbitraria con (digamos) media $(\mu_x,\mu_y),$ desviaciones típicas positivas $(\sigma_x, \sigma_y),$ y correlación $\rho.$ Como no se trata de nada nuevo, sino simplemente de cambiar las unidades de medida de las dos variables, me limitaré a exponer los resultados.
La densidad del máximo $Z=\max(X,Y)$ viene dada por la suma de dos términos. La primera (en la que $X \ge Y$ ) contribuye
$$\phi(z;\ \mu_x, \sigma_x)\ \Phi\left(z;\ h(z; \mu, \sigma, \rho),\ \sigma_y \sqrt{1-\rho^2}\right)$$
donde $\phi$ es la densidad Normal (con media y desviación típica dadas), $\Phi$ es la función de distribución Normal (con los mismos significados de los parámetros), y
$$h(z; \mu, \sigma, \rho) = E[Y\mid X=z] = \mu_y + \frac{\sigma_y}{\sigma_x}\, \rho\, (z - \mu_x)$$
es la expectativa condicional de la segunda variable dada la primera igual a $z.$
El segundo término se obtiene cambiando los papeles de las variables (de modo que $Y\ge X$ ): aunque $z$ (que representa su máximo) permanece igual en la fórmula y $\rho$ tampoco cambia, $\mu$ y $\sigma$ se invierten.
Aquí, para ser totalmente explícitos, hay una R
aplicación de la función de densidad de $Z$ basándose directamente en esta descripción.
#
# Density of the maximum of a nondegenerate bivariate Normal.
#
dnorm2max <- function(z, mu=c(0,0), sigma=c(1,1), rho=0) {
# The conditional expectation at `z`.
h <- function(z, m, s) m[2] + s[2] / s[1] * rho * (z - m[1])
# The contribution of one (infinitesimal) strip.
f <- function(z, m, s) dnorm(z, m[1], s[1]) * pnorm(z, h(z, m, s), s[2] * sqrt(1 - rho^2))
# Return the total from the vertical and horizontal strips.
f(z, mu, sigma) + f(z, rev(mu), rev(sigma))
}
Como un cheque, esta figura muestra un histograma de $Z$ de una simulación de longitud 100.000, en la que he superpuesto el gráfico de dnorm2max
. El acuerdo es excelente.
Este ejemplo concreto ayuda a comprenderlo mejor y ayuda a informar nuestra intuición. Porque $b\approx -1,$ esto es un fuerte anticorrelacionado serie: sus valores tienden a oscilar entre cifras por encima y por debajo de la media.
La varianza de la primera diferencia $x_{t+1}-x_t,$ igual a $50,$ es relativamente grande. La diferencia de dos pasos $x_{t+2}-x_t$ sólo tiene una varianza de $2,$ lo que indica que cada primera diferencia se invierte casi por completo en el paso siguiente.
La correlación entre estas dos variables -las diferencias de un paso y de dos pasos- es baja. $0.1,$ sin embargo. En consecuencia, es muy probable que el máximo de las dos diferencias sea la primera cuando es positiva y, en caso contrario, es probable que sea la segunda cuando la primera es negativa. Por lo tanto, el resultado se parece un poco a una mezcla igual de una variable Normal con varianza $2$ (verde en la figura siguiente) y una variable seminormal con varianza $50$ (azul).
Por fin, una forma atractiva de estimar la distribución de $Z$ (o cualquiera de sus propiedades, como los percentiles) sería estimar el parámetro $b$ de la serie $(x_t).$ (Siendo realistas, es probable que tenga que calcular $\sigma^2$ también). Puede propagar los errores estándar utilizando el método Delta. Cuando sus estimaciones se basan en la optimización de una probabilidad, este método está rigurosamente justificado.
Con una realización suficientemente larga de este modelo de series temporales se podría, por supuesto, estimar la distribución de $Z$ a partir de su distribución empírica. Puede que merezca la pena comprobar sus estimaciones.