Simulation eines Gaußschen Prozesses (Ornstein Uhlenbeck) mit einer exponentiell abfallenden Kovarianzfunktion

Ich versuche, viele Draws (dh Realisierungen) eines Gaußschen Prozesses , mit dem Mittelwert 0 und der Kovarianzfunktion zu erzeugen . $e_i(t)$ $1\leq t \leq T$ $\gamma(s,t)=\exp(-|t-s|)$

Gibt es einen effizienten Weg, um die Quadratwurzel einer Kovarianzmatrix nicht zu berechnen? Kann jemand alternativ ein Paket empfehlen , um dies zu tun? $T \times T$ R

r probability simulation stochastic-processes gaussian-process user603
quelle

Es ist ein stationärer Prozess (ähnelt einer einfachen Version eines OU-Prozesses). Wird es gleichmäßig abgetastet?

Kardinal

Das R-Paket mvtnormhat rmvnorm(n, mean, sigma)wo sigmaist die Kovarianzmatrix; Sie müssten jedoch die Kovarianzmatrix für Ihre abgetasteten / ausgewählten selbst erstellen.

t

$t$

Jbowman

@jb Vermutlich ist riesig, sonst würde die OP nicht fragen die Matrix Zersetzung zu vermeiden (was implizit in ist ).

T

$T$ rmvnorm

whuber

@cardinal Ich stimme zu, dies ist ein Ornstein-Uhlenbeck-Gauß-Prozess. (Es wäre großartig, wenn das Schlüsselwort "Ornstein Uhlenbeck" in die Frage und / oder den Titel umgewandelt werden könnte. Diese Frage würde umso mehr Verkehr erhalten, je mehr Verkehr sie verdient.)

redmoskito

Antworten:

Ja. Es gibt einen sehr effizienten (linearen Zeit-) Algorithmus, und die Intuition dafür stammt direkt aus dem einheitlich abgetasteten Fall.

Angenommen , wir haben eine Partition von , so dass . $[0,T]$ $0=t_0 < t_1 < t_2 < \cdots < t_n = T$

Fall mit einheitlicher Stichprobe

In diesem Fall haben wir wobei . Let den Wert des diskret abgetasteten Prozess zum Zeitpunkt bezeichnen . $t_i = i \Delta$ $\Delta = T/n$ $X_i := X(t_i)$ $t_i$

Es ist leicht zu erkennen, dass einen AR (1) -Prozess mit der Korrelation . Daher können wir einen Beispielpfad für die Partition wie folgt erzeugen: wobei sind iid und . $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$

Allgemeiner Fall

Wir könnten uns dann vorstellen, dass dies für eine allgemeine Partition möglich sein könnte. Insbesondere sei und . Wir haben das und so könnten wir vermuten, dass $\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} \>.$

In der Tat ist und so haben wir zumindest die Korrelation mit dem benachbarten Term korrekt. $\mathbb E X_{i+1} X_i = \rho_i$

Das Ergebnis folgt nun durch Teleskopieren über die Turmeigenschaft der bedingten Erwartung. Nämlich und die Produktteleskope in der auf folgende Weise

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

Dies beweist das Ergebnis. Daher kann der Prozess auf einer beliebigen Partition aus einer Folge von iid Zufallsvariablen in -Zeit erzeugt werden, wobei die Größe der Partition ist. $\mathcal N(0,1)$ $O(n)$ $n$

NB : Dies ist eine exakte Abtasttechnik, da sie eine abgetastete Version des gewünschten Prozesses mit den genau korrekten endlichdimensionalen Verteilungen liefert . Dies steht im Gegensatz zu Euler- (und anderen) Diskretisierungsschemata für allgemeinere SDEs, bei denen aufgrund der Annäherung über die Diskretisierung eine Verzerrung auftritt.

Kardinal
quelle

Nur noch ein paar Bemerkungen. (1) Um eine gute Vorstellung davon zu bekommen, wie der kontinuierliche Zeitprozess aussieht, müssen und so gewählt werden, dass klein ist, sagen wir weniger als . (2) Die inverse Kovarianzmatrix ( Präzisionsmatrix ) für den Zeitreihenvektor ist dreidiagonal, ebenso wie seine Cholesky-Wurzel.

n

$n$

T

$T$

Δ

$\Delta$

0.1

$0.1$

Yves

@ Yves: Danke für deine Kommentare. Um es klar auszudrücken, bietet das von mir beschriebene Verfahren eine genaue Realisierung des zeitkontinuierlichen Prozesses, der auf der entsprechenden Partition abgetastet wurde. Insbesondere gibt es keinen Diskretisierungsfehler wie bei einer typischen Euler-Schema-Annäherung an allgemeinere SDEs. Der inverse Cholesky hat, wie die Konstruktion in der Antwort zeigt, nur in der Diagonale und in der unteren Off-Diagonale Begriffe ungleich Null, ist also etwas einfacher als die Tridiagonale.

Kardinal

Gute Antwort! Verallgemeinert sich dies auf den allgemeinen OU-Prozess mit beliebiger Skala, ? Es scheint so.

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

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

Redmoskito

Berechnen Sie die zerlegte Kovarianzmatrix durch unvollständige Cholesky-Zerlegung oder eine andere Matrixzerlegungstechnik. Die zerlegte Matrix sollte TxM sein, wobei M nur ein Bruchteil von T ist.

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

Steven
quelle

Können Sie hier eine explizite Form der Cholesky-Zerlegung angeben? Ich denke, dass die Antwort des Kardinals genau das erreicht, wenn Sie darüber nachdenken, indem Sie als Funktion der Geschichte ausdrücken .

X_{i}

$X_i$

StasK

Der Algorithmus ist etwas zu lang, um ihn zusammenzufassen. Eine hervorragende Beschreibung finden Sie hier: Kernel ICA , Seite 20. Beachten Sie, dass dieser Algorithmus unvollständig ist , dh nicht die gesamte Zerlegung berechnet, sondern eine Annäherung (daher ist er viel schneller). Ich habe Code für diesen Algorithmus in der KMBOX-Toolbox veröffentlicht. Sie können ihn hier herunterladen: km_kernel_icd .

Steven