Simulação de um processo gaussiano (Ornstein Uhlenbeck) com uma função de covariância exponencialmente decadente

8

Eu estou tentando gerar muitos draws (ou seja, realizações) de um processo de Gauss ei(t) , 1tT com média 0 e função covariância γ(s,t)=exp(|ts|) .

Existe uma maneira eficiente de fazer isso que não envolveria calcular a raiz quadrada de um T×T matriz de covariância? Como alternativa, alguém pode recomendar um Rpacote para fazer isso?

user603
fonte
2
É um processo estacionário (parece próximo a uma versão simples de um processo de OU). É amostrada uniformemente?
cardeal
O pacote R mvtnormpossui rmvnorm(n, mean, sigma)onde sigmaestá a matriz de covariância; você teria que construir a matriz de covariância para o seu amostrados / selecionada t é mesmo, embora.
jbowman
2
@jb Presumivelmente T é enorme, caso contrário , o OP não pediria para evitar a decomposição da matriz (que está implícita rmvnorm).
whuber
1
@ cardinal Concordo, este é um processo Gaussiano de Ornstein-Uhlenbeck. (Seria ótimo se a palavra-chave "Ornstein Uhlenbeck" pode ser editado na questão e / ou título Ficaria esta questão, o mais tráfego que merece.)
redmoskito

Respostas:

11

Sim. Existe um algoritmo muito eficiente (tempo linear), e a intuição para isso vem diretamente do caso de amostra uniforme.

Suponhamos que temos uma partição de de modo a que 0 = t 0 < t 1 < t 2 < < t n = T .[0,T]0=t0<t1<t2<<tn=T

Caso de amostra uniforme

Neste caso, temos onde Δ = T / n . Seja X i : = X ( t i ) denota o valor do processo amostrado discretamente no momento t i .ti=iΔΔ=T/nXi:=X(ti)ti

É fácil de ver que o formar um processo AR (1) com correlação ρ = exp ( - Δ ) . Portanto, podemos gerar um caminho de amostra { X t } para a partição da seguinte forma X i + 1 = ρ X i + Xiρ=exp(Δ){Xt} Em que Z i são iid N ( 0 , 1 ) e X 0 = Z 0 .

Xi+1=ρXi+1ρ2Zi+1,
ZiN(0,1)X0=Z0

Caso Geral

Podemos então imaginar que seria possível fazer isso para uma partição geral . Em particular, seja e ρ i = exp ( - Δ i ) . Temos que γ ( t i , t i + 1 ) = ρ iΔi=ti+1tiρi=exp(Δi) e , portanto, podemos supor que X i + 1 = ρ i X i +

γ(ti,ti+1)=ρi,
Xi+1=ρiXi+1ρi2Zi+1.

EXi+1Xi=ρi

EXiXi=E(E(XiXiXi1))=ρi1EXi1Xi==k=1ρik,
k=1ρik=exp(k=1Δik)=exp(titi)=γ(ti,ti).

N(0,1)O(n)n

NB : Esta é uma técnica de amostragem exata , pois fornece uma versão amostrada do processo desejado com as distribuições de dimensões finitas exatamente corretas . Isso contrasta com os esquemas de discretização de Euler (e outros) para SDEs mais gerais, que sofrem um viés devido à aproximação via discretização.

cardeal
fonte
nTΔ0.1
@ Yves: Obrigado por seus comentários. Para ser claro, o procedimento que descrevi fornece uma realização exata do processo de tempo contínuo amostrado na partição correspondente; em particular, não erro de discretização, como ocorre na aproximação típica do esquema de Euler a SDEs mais gerais. O Cholesky inverso, como mostra a construção na resposta, tem termos diferentes de zero apenas na diagonal e na diagonal inferior, portanto é um pouco mais simples que o tridiagonal.
cardeal
γ(ti,tj)=exp(α|titj|)
1

Calcule a matriz de covariância decomposta por decomposição incompleta de Cholesky ou qualquer outra técnica de decomposição de matriz. A matriz decomposta deve ser TxM, onde M é apenas uma fração de T.

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

Steven
fonte
2
Xi
1
O algoritmo é um pouco longo para resumir. Você pode encontrar uma excelente descrição aqui: Kernel ICA , página 20. Observe que esse algoritmo está incompleto , o que significa que não calcula toda a decomposição, mas uma aproximação (portanto, é muito mais rápido). Publiquei o código desse algoritmo na caixa de ferramentas do KMBOX, você pode baixá-lo aqui: km_kernel_icd .
6303 Steven