Simuler un processus gaussien (Ornstein Uhlenbeck) avec une fonction de covariance à décroissance exponentielle

J'essaie de générer de nombreux tirages (c'est-à-dire des réalisations) d'un processus gaussien , avec la moyenne 0 et la fonction de covariance . $e_i(t)$ $1\leq t \leq T$ $\gamma(s,t)=\exp(-|t-s|)$

Existe-t-il un moyen efficace de le faire qui n'impliquerait pas de calculer la racine carrée d'une matrice de covariance ? Sinon, quelqu'un peut-il recommander un package pour ce faire? $T \times T$ R

r probability simulation stochastic-processes gaussian-process user603
la source

C'est un processus stationnaire (ressemble à une version simple d'un processus OU). Est-il échantillonné uniformément?

cardinal

Le package R mvtnorma rmvnorm(n, mean, sigma)où sigmaest la matrice de covariance; vous devez cependant construire vous-même la matrice de covariance pour vos

échantillonnés / sélectionnés

t

$t$ .

jbowman

@jb Vraisemblablement est énorme, sinon l'OP ne demanderait pas à éviter la décomposition de la matrice (qui est implicite dans ).

T

$T$ rmvnorm

whuber

@cardinal Je suis d'accord, c'est un processus gaussien Ornstein-Uhlenbeck. (Ce serait génial si le mot-clé "Ornstein Uhlenbeck" pouvait être modifié dans la question et / ou le titre. Il obtiendrait cette question plus de trafic qu'elle mérite)

redmoskito

Réponses:

Oui. Il existe un algorithme très efficace (temps linéaire), et son intuition vient directement du cas échantillonné uniformément.

Supposons que nous ayons une partition de tel que . $[0,T]$ $0=t_0 < t_1 < t_2 < \cdots < t_n = T$

Cas échantillonné uniformément

Dans ce cas, nous avons où . Soit la valeur du processus échantillonné discrètement au temps . $t_i = i \Delta$ $\Delta = T/n$ $X_i := X(t_i)$ $t_i$

Il est facile de voir que les forment un processus AR (1) avec corrélation . Par conséquent, nous pouvons générer un exemple de chemin pour la partition comme suit où sont iid et . $X_i$ $\rho = \exp(-\Delta)$ $\{X_t\}$

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

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

Z_{i}

$Z_i$

N (0, 1)

$\mathcal N(0,1)$

X_{0} = Z_{0}

$X_0 = Z_0$

Cas général

On pourrait alors imaginer qu'il pourrait être possible de le faire pour une partition générale . En particulier, soit et . Nous avons cela et nous pouvons donc deviner que $\Delta_i = t_{i+1} - t_i$ $\rho_i = \exp(-\Delta_i)$

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

$\gamma(t_i,t_{i+1}) = \rho_i \>,$

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 effet, et donc nous avons au moins la corrélation avec le terme voisin correcte. $\mathbb E X_{i+1} X_i = \rho_i$

Le résultat suit maintenant en se télescopant via la propriété de la tour de l'attente conditionnelle. A savoir, et les télescopes de produit dans le de la manière suivante

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} \>,$

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

Cela prouve le résultat. Par conséquent, le processus peut être généré sur une partition arbitraire à partir d'une séquence de variables aléatoires iid en temps où est la taille de la partition. $\mathcal N(0,1)$ $O(n)$ $n$

NB : Il s'agit d'une technique d'échantillonnage exacte en ce qu'elle fournit une version échantillonnée du processus souhaité avec les distributions de dimensions finies exactement correctes . Cela contraste avec les schémas de discrétisation d'Euler (et d'autres) pour les SDE plus généraux, qui encourent un biais en raison de l'approximation via la discrétisation.

cardinal
la source

Encore quelques remarques. (1) Pour avoir une bonne idée de ce à quoi ressemble le processus de temps continu, et doivent être choisis de sorte que soit petit, disons inférieur à . (2) La matrice de covariance inverse ( précision ) pour le vecteur série temporelle est tri-diagonale, tout comme sa racine de Cholesky.

n

$n$

T

$T$

Δ

$\Delta$

0.1

$0.1$

Yves

@Yves: Merci pour vos commentaires. Pour être clair, la procédure que j'ai décrite donne une réalisation exacte du processus en temps continu échantillonné sur la partition correspondante; en particulier, il n'y a pas d' erreur de discrétisation comme dans l'approximation du schéma d'Euler typique des SDE plus généraux. L'inverse de Cholesky, comme le montre la construction de la réponse, n'a des termes non nuls que sur la diagonale et la diagonale inférieure, c'est donc un peu plus simple que tridiagonal.

cardinal

Très bonne réponse! Est-ce que cela se généralise au processus OU général avec une échelle arbitraire, ? Il semble que ce soit possible.

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

$\gamma(t_i, t_j) = \exp(\alpha \; |t_i - t_j|)$

redmoskito

Calculez la matrice de covariance décomposée par décomposition Cholesky incomplète ou toute autre technique de décomposition matricielle. La matrice décomposée doit être TxM, où M n'est qu'une fraction de T.

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

Steven
la source

Pouvez-vous donner ici une forme explicite de la décomposition de Cholesky? Je pense que la réponse du cardinal y parvient, si vous y réfléchissez, en exprimant en fonction de l'histoire.

X_{i}

$X_i$

StasK

L'algorithme est un peu trop long pour résumer. Vous pouvez trouver une excellente description ici: Kernel ICA , page 20. Notez que cet algorithme est incomplet , ce qui signifie qu'il ne calcule pas la décomposition entière mais plutôt une approximation (il est donc beaucoup plus rapide). J'ai publié le code de cet algorithme dans la boîte à outils KMBOX, vous pouvez le télécharger ici: km_kernel_icd .

Steven