【问题标题】:Simple Gaussian Process Simulation on RR上的简单高斯过程模拟
【发布时间】:2020-10-31 16:50:50
【问题描述】:

如何模拟高斯过程X(t), t = 1, . . . , 200,均值函数m(t) = 0 和 协方差函数r(h) = Cov(t, t + h) = exp(-|h|)。我知道这个过程有时被称为 Ornstein-Uhlenbeck 过程,但如何绘制模拟过程。

感谢期待

【问题讨论】:

  • 到目前为止您尝试过什么?我相信维基百科页面上显示的连续时间方程的离散化应该会有所帮助:X(t+1) = h*X(t) + rnorm(1) 应该可以工作
  • 您知道具有该协方差结构的过程的公式吗?
  • 是的,它的:r(h) = exp(-|h|)
  • 具有该协方差结构的过程的公式是 r(h) = exp(-|h|) @ Ricardo Semião e Castro

标签: r stochastic stochastic-process


【解决方案1】:

按照Wikipedia 的定义,Ornstein–Uhlenbeck 过程由以下随机微分方程定义: 其中 是一个维纳过程,它的一个属性是它具有高斯增量,即

上述方程可以用以下方式离散化: 在哪里 由于维纳过程的高斯增量属性,我们有。这意味着可以使用sqrt(dt)*rnorm(1)生成增量值

我在 R 中编写了以下函数,它采用时间向量、过程平均值、标准差和 theta 值。

simulate <- function(t, mean=0, std=1, x0=mean, theta=1, number.of.points=length(t)){
  # calculate time differences
  dt <- diff(t)
  X <- vector("numeric", length=number.of.points)
  X[1] <- x0
  for(i in 1:(number.of.points-1)){
    X[i+1] <- X[i] + theta * (mean-X[i])*dt[i] + std * sqrt(dt[i])* rnorm(1)
  }
  data.frame(x=t, y=X)
}
simulate(t=1:200) %>%  ggplot(aes(x,y)) + geom_line()

另一个使用purrr的实现

simulate <- function(t, mean=0, sd=1, theta=1, number.of.points=length(t)){
  stopifnot(!missing(t) | !missing(number.of.points))
  if(missing(t)){
    t <- 1:number.of.points
  } 

  unlist(purrr::accumulate2(vector("numeric", length=number.of.points-1), diff(t), function(x, o, y) {
    x + theta*(mean - x)* y + sqrt(y)*rnorm(1)
  }, .init=x0), use.names=F) -> X
  
  data.frame(x=t, y=X)
}
simulate(number.of.points=200) %>%  ggplot(aes(x,y)) + geom_line()

【讨论】:

  • 非常感谢。 @Abdessabour Mtk,我必须在这里管理一种稍微不同的方法,因为如何从基本构建块中做到这一点,例如通过计算 X = (X1, X2 ... X200) 的均值和协方差矩阵,然后随机绘制使用函数“rmnorm”来自多元正态分布的向量?那怎么办?有什么手续吗?
  • @Ticktock 看看这个answer 可能会很有趣,如果你使用答案,我认为你需要运行多个模拟才能计算均值和协方差矩阵使用rowMeans(sims) 标记将为您提供方法,cov(t(sims)) 将为您提供协方差矩阵
  • @Abdessabour Mtk 我仍然无法从基本构建块中生成结果,请您在计算 X = (X1, X2 ... X200) 的均值和协方差矩阵时提供帮助然后使用函数“rmnorm”从多元正态分布中绘制一个随机向量?协方差函数也是r(h)=exp(-|h|)。非常感谢您的期待,先生!
【解决方案2】:

从这里使用函数:https://quant.stackexchange.com/questions/1260/r-code-for-ornstein-uhlenbeck-process

ornstein_uhlenbeck <- function(T,n,nu,lambda,sigma,x0){
  dw  <- rnorm(n, 0, sqrt(T/n))
  dt  <- T/n
  x <- c(x0)
  for (i in 2:(n+1)) {
    x[i]  <-  x[i-1] + lambda*(nu-x[i-1])*dt + sigma*dw[i-1]
  }
  return(x);
}

test <- ornstein_uhlenbeck(200, 200, 0, 0.8, 1, 0)
plot(x = seq_along(test), y = test, type = 'l')

(请注意,它只为您提供链接中问题答案之一所述的近似分布。)

我假设 T = 200, n = 200, nu = 0(如你所述),均值回归参数为 0.8,sigma 为 1,过程从 0 开始。

【讨论】:

  • 非常感谢大家。 @tester,我必须在这里管理一种稍微不同的方法,即如何从基本构建块中做到这一点,例如通过计算 X = (X1, X2 ... X200) 的均值和协方差矩阵,然后绘制随机向量从使用函数“rmnorm”的多元正态分布?
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2013-12-16
  • 1970-01-01
  • 2015-07-12
  • 1970-01-01
  • 2014-04-17
  • 1970-01-01
  • 2015-01-02
相关资源
最近更新 更多