Simulación de un proceso gaussiano (Ornstein Uhlenbeck) con una función de covarianza exponencialmente decadente

Estoy tratando de generar muchos sorteos (es decir, realizaciones) de un proceso gaussiano $e_i(t)$ , $1\leq t \leq T$ con media 0 y función de covarianza $\gamma(s,t)=\exp(-|t-s|)$ .

¿Hay una manera eficiente de hacer esto que no implique el cálculo de la raíz cuadrada de un $T \times T$ ¿Matriz de covarianza? Alternativamente, ¿alguien puede recomendar un Rpaquete para hacer esto?

r probability simulation stochastic-processes gaussian-process usuario603
fuente

Es un proceso estacionario (se parece a una versión simple de un proceso OU). ¿Se muestrea de manera uniforme?

cardenal

El paquete R mvtnormtiene rmvnorm(n, mean, sigma)dónde sigmaestá la matriz de covarianza; tendrías que construir la matriz de covarianza para tu muestra / seleccionado

t

$t$ es usted mismo, sin embargo.

jbowman

@jb Presumiblemente

T

$T$ es enorme, de lo contrario , el OP no estaría pidiendo evitar la descomposición de la matriz (que está implícita en rmvnorm).

whuber

@cardinal Estoy de acuerdo, este es un proceso gaussiano de Ornstein-Uhlenbeck. (Sería genial si la palabra clave "Ornstein Uhlenbeck" se pudiera editar en la pregunta y / o título.

Obtendría

Respuestas:

Si. Existe un algoritmo muy eficiente (tiempo lineal), y su intuición proviene directamente del caso de muestra uniforme.

Supongamos que tenemos una partición de $[0,T]$ tal que $0=t_0 < t_1 < t_2 < \cdots < t_n = T$ .

Estuche uniformemente muestreado

En este caso tenemos $t_i = i \Delta$ dónde $\Delta = T/n$ . Dejar $X_i := X(t_i)$ denotar el valor del proceso discretamente muestreado en el momento $t_i$ .

Es fácil ver que el $X_i$ formar un proceso AR (1) con correlación $\rho = \exp(-\Delta)$ . Por lo tanto, podemos generar una ruta de muestra $\{X_t\}$ para la partición de la siguiente manera

X_{i + 1} = ρ X_{i} + \sqrt{1 - ρ^{2}} Z_{i + 1},

$X_{i+1} = \rho X_i + \sqrt{1-\rho^2} Z_{i+1} \>,$ dónde

Z_{i}

$Z_i$ son iid

N (0, 1)

$\mathcal N(0,1)$ y

X_{0} = Z_{0}

$X_0 = Z_0$ .

Caso general

Entonces podríamos imaginar que podría ser posible hacer esto para una partición general . En particular, dejemos $\Delta_i = t_{i+1} - t_i$ y $\rho_i = \exp(-\Delta_i)$ . Tenemos eso

γ (t_{i}, t_{i + 1}) = ρ_{i},

$\gamma(t_i,t_{i+1}) = \rho_i \>,$ y entonces podríamos adivinar que

X_{i + 1} = ρ_{i} X_{i} + \sqrt{1 - ρ_{i}^{2}} Z_{i + 1} .

$X_{i+1} = \rho_i X_i + \sqrt{1-\rho_i^2} Z_{i+1} \>.$

En efecto, $\mathbb E X_{i+1} X_i = \rho_i$ y así al menos tenemos la correlación correcta con el término vecino.

El resultado ahora se obtiene mediante la telescopía a través de la propiedad de la torre de expectativa condicional. A saber,

E X_{i} X_{i - ℓ} = E (E (X_{i} X_{i - ℓ} ∣ X_{i - 1})) = ρ_{i - 1} E X_{i - 1} X_{i - ℓ} = \dots = \prod_{k = 1}^{ℓ} ρ_{i - k},

$\newcommand{\e}{\mathbb E} \e X_i X_{i-\ell} = \e( \e(X_i X_{i-\ell} \mid X_{i-1} )) = \rho_{i-1} \mathbb E X_{i-1} X_{i-\ell} = \cdots = \prod_{k=1}^\ell \rho_{i-k} \>,$ y los telescopios del producto de la siguiente manera

\prod_{k = 1}^{ℓ} ρ_{i - k} = \exp (- \sum_{k = 1}^{ℓ} Δ_{i - k}) = \exp (t_{i - ℓ} - t_{i}) = γ (t_{i - ℓ}, t_{i}) .

$\prod_{k=1}^\ell \rho_{i-k} = \exp\Big(-\sum_{k=1}^\ell \Delta_{i-k}\Big) = \exp(t_{i-\ell} - t_i) = \gamma(t_{i-\ell},t_i) \>.$

Esto prueba el resultado. Por lo tanto, el proceso se puede generar en una partición arbitraria a partir de una secuencia de iid $\mathcal N(0,1)$ variables aleatorias en $O(n)$ tiempo donde $n$ es el tamaño de la partición

NB : Esta es una técnica de muestreo exacta , ya que proporciona una versión muestreada del proceso deseado con las distribuciones finitas dimensionales correctas . Esto contrasta con los esquemas de discretización de Euler (y otros) para las SDE más generales, que incurren en un sesgo debido a la aproximación por discretización.

cardenal
fuente

Solo unas pocas observaciones más. (1) Para tener una buena idea de cómo es el proceso de tiempo continuo,

n

$n$ y

T

$T$ debe ser elegido para que

Δ

$\Delta$ es pequeño, digamos menos de

0.1

$0.1$ . (2) La matriz de covarianza inversa ( precisión ) para el vector de series temporales es tri-diagonal, al igual que su raíz de Cholesky.

Yves

@Yves: Gracias por tus comentarios. Para ser claros, el procedimiento que describí da una realización exacta del proceso de tiempo continuo muestreado en la partición correspondiente; en particular, no hay error de discretización como existe en la aproximación típica del esquema de Euler a las SDE más generales. El Cholesky inverso, como lo muestra la construcción en la respuesta, tiene términos distintos de cero solo en diagonal y fuera de diagonal inferior, por lo que es un poco más simple que tridiagonal.

cardenal

¡Gran respuesta! ¿Esto se generaliza al proceso general de OU con escala arbitraria,

γ (t_{i}, t_{j}) = \exp (α | t_{i} - t_{j} |)

$\gamma(t_i, t_j) = \exp(\alpha \; |t_i - t_j|)$ ? Parece que podría ser.

redmoskito

Calcule la matriz de covarianza descompuesta por descomposición de Cholesky incompleta o cualquier otra técnica de descomposición de matriz. La matriz descompuesta debe ser TxM, donde M es solo una fracción de T.

http://en.wikipedia.org/wiki/Incomplete_Cholesky_factorization

Steven
fuente

¿Puedes dar una forma explícita de la descomposición de Cholesky aquí? Creo que la respuesta del cardenal logra exactamente eso, si lo piensas, expresando

X_{i}

$X_i$ en función de la historia.

StasK

El algoritmo es demasiado largo para resumirlo. Puede encontrar una excelente descripción aquí: Kernel ICA , página 20. Tenga en cuenta que este algoritmo está incompleto , lo que significa que no calcula la descomposición completa sino más bien una aproximación (por lo tanto, es mucho más rápido). He publicado el código para este algoritmo en la caja de herramientas KMBOX, puede descargarlo aquí: km_kernel_icd .

Steven