Mensimulasikan proses Gaussian (Ornstein Uhlenbeck) dengan fungsi kovarians yang membusuk secara eksponensial

Saya mencoba untuk menghasilkan banyak undian (yaitu, realisasi) dari proses Gaussian $e_i(t)$ , $1\leq t \leq T$ dengan rata-rata 0 dan fungsi kovarian $\gamma(s,t)=\exp(-|t-s|)$ .

Apakah ada cara yang efisien untuk melakukan ini yang tidak akan melibatkan perhitungan akar kuadrat dari a $T \times T$ matriks kovarians? Atau adakah yang bisa merekomendasikan Rpaket untuk melakukan ini?

r probability simulation stochastic-processes gaussian-process pengguna603
sumber

Ini adalah proses stasioner (terlihat dekat dengan versi sederhana dari proses OU). Apakah sampelnya seragam?

kardinal

Paket R mvtnormmemiliki di rmvnorm(n, mean, sigma)mana sigmamatriks kovarians; Anda harus membuat matriks kovarians untuk sampel / yang Anda pilih

t

$t$ itu sendiri, meskipun.

Jbowman

@ jb Agaknya

T

$T$ sangat besar, jika OP tidak akan meminta untuk menghindari dekomposisi matriks (yang tersirat dalam rmvnorm).

whuber

@ cardinal Saya setuju, ini adalah proses Gaussian Ornstein-Uhlenbeck. (Alangkah baiknya jika kata kunci "Ornstein Uhlenbeck" dapat diedit menjadi pertanyaan dan / atau judul. Akan mendapatkan pertanyaan ini semakin banyak lalu lintas yang layak)

redmoskito

Jawaban:

Iya. Ada algoritma (waktu linear) yang sangat efisien, dan intuisi untuknya datang langsung dari kasing yang seragam.

Misalkan kita memiliki partisi $[0,T]$ seperti yang $0=t_0 < t_1 < t_2 < \cdots < t_n = T$ .

Kasing sampel yang seragam

Dalam hal ini yang kita miliki $t_i = i \Delta$ dimana $\Delta = T/n$ . Membiarkan $X_i := X(t_i)$ menunjukkan nilai dari proses sampel yang diambil pada saat itu $t_i$ .

Sangat mudah untuk melihat bahwa $X_i$ membentuk proses AR (1) dengan korelasi $\rho = \exp(-\Delta)$ . Karenanya, kita dapat menghasilkan jalur sampel $\{X_t\}$ untuk partisi sebagai berikut

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

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

Z_{i}

$Z_i$ apakah iid

N (0, 1)

$\mathcal N(0,1)$ dan

X_{0} = Z_{0}

$X_0 = Z_0$ .

Kasus umum

Kita kemudian dapat membayangkan bahwa ini dapat dilakukan untuk partisi umum . Secara khusus, mari $\Delta_i = t_{i+1} - t_i$ dan $\rho_i = \exp(-\Delta_i)$ . Kami memilikinya

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

$\gamma(t_i,t_{i+1}) = \rho_i \>,$ dan jadi kita mungkin menebaknya

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

Memang, $\mathbb E X_{i+1} X_i = \rho_i$ dan setidaknya kita memiliki korelasi dengan istilah tetangga yang benar.

Hasilnya sekarang mengikuti dengan telescoping melalui properti menara harapan bersyarat. Yaitu, dan teleskop produk di dengan cara berikut

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

Ini membuktikan hasilnya. Oleh karena itu proses dapat dihasilkan pada partisi arbitrer dari urutan variabel acak iid dalam waktu di mana adalah ukuran partisi. $\mathcal N(0,1)$ $O(n)$ $n$

NB : Ini adalah teknik pengambilan sampel yang tepat karena memberikan versi sampel dari proses yang diinginkan dengan distribusi dimensi-terbatas yang benar - benar tepat . Ini berbeda dengan skema diskritisasi Euler (dan lainnya) untuk SDE yang lebih umum, yang menimbulkan bias karena perkiraan melalui diskretisasi.

kardinal
sumber

Hanya beberapa komentar lagi. (1) Untuk mendapatkan gambaran yang baik tentang seperti apa proses waktu kontinu, dan harus dipilih sehingga kecil, katakan kurang dari . (2) Matriks kovarians terbalik ( presisi ) untuk vektor deret waktu adalah tri-diagonal, seperti halnya akar Cholesky-nya.

n

$n$

T

$T$

Δ

$\Delta$

0.1

$0.1$

Yves

@Yves: Terima kasih atas komentar Anda. Agar lebih jelas, prosedur yang telah saya uraikan memberikan realisasi yang tepat dari proses waktu kontinu yang disampel pada partisi terkait; khususnya, tidak ada kesalahan diskritisasi seperti ada dalam pendekatan skema-Euler yang khas untuk SDE yang lebih umum. Pembalikan Cholesky, seperti yang ditunjukkan oleh konstruksi dalam jawaban memiliki istilah bukan nol hanya pada diagonal dan lebih rendah di-diagonal, sehingga sedikit lebih sederhana daripada tridiagonal.

kardinal

Jawaban bagus! Apakah ini digeneralisasi ke proses OU umum dengan skala arbitrer, ? Sepertinya itu mungkin.

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

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

redmoskito

Hitung matriks kovarians terurai dengan dekomposisi Cholesky tidak lengkap atau teknik dekomposisi matriks lainnya. Matriks yang diuraikan harus TxM, di mana M hanya sebagian kecil dari T.

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

Steven
sumber

Bisakah Anda memberikan bentuk eksplisit dari dekomposisi Cholesky di sini? Saya pikir jawaban kardinal mencapai hal itu, jika Anda memikirkannya, dengan menyatakan sebagai fungsi sejarah.

X_{i}

$X_i$

Tugas

Algoritma ini agak terlalu panjang untuk diringkas. Anda dapat menemukan deskripsi yang sangat baik di sini: Kernel ICA , halaman 20. Perhatikan bahwa algoritma ini tidak lengkap , artinya tidak menghitung keseluruhan dekomposisi tetapi lebih merupakan perkiraan (karenanya jauh lebih cepat). Saya telah menerbitkan kode untuk algoritma ini di kotak alat KMBOX, Anda dapat mengunduhnya di sini: km_kernel_icd .

Steven