两个布朗运动之间相关性的蒙特卡罗模拟（连续随机游走）答案

【问题标题】：Monte Carlo simulation of correlation between two Brownian motion (continuous random walk)两个布朗运动之间相关性的蒙特卡罗模拟（连续随机游走）
【发布时间】：2017-02-18 00:58:37
【问题描述】：

y <- cumsum(rnorm(100,0,1)) # random normal, with small (1.0) drift.
y.ts <- ts(y)
x <- cumsum(rnorm(100,0,1))
x
x.ts <- ts(x)
ts.plot(y.ts,ty= "l", x.ts) # plot the two random walks 


Regression.Q1 = lm(y~x) ; summary(lm2) 
summary(Regression.Q1)

t.test1 <- (summary(Regression.Q1)$coef[2,3]) # T-test computation 


y[t] = y[t-1] + epsilon[t]
epsilon[t] ~ N(0,1)
set.seed(1)
t=1000
epsilon=sample(c(-1,1), t, replace = 1) # Generate k random walks across time {0, 1, ... , T}


N=T=1e3
y=t(apply(matrix(sample(c(-1,1),N*T,rep=TRUE),ncol=T),1,cumsum))
y[1]<-0
for (i in 2:t) {
  y[i]<-y[i-1]+epsilon[i]
}

我需要：

重复该过程 1000 次（蒙特卡洛模拟），即围绕前一个程序构建一个循环，每次保存 t 统计量。您将有 1;000 个 t-tests 序列：S = (t-test1, t-test2, ... , t-test1000)。计算 1,000 t 检验的绝对值 > 1.96 的次数，即 5% 显着性水平的临界值。如果系列是 I(0)，你会发现大约 5%。这里不是这种情况（虚假回归）。

我需要添加什么来保存各自的系数？

【问题讨论】：

标签： r loops simulation montecarlo

【解决方案1】：

您发布的与y[t] = y[t-1] + epsilon[t] 相关的代码不是真正的工作代码，但我可以看到您正在尝试存储所有 1000 * 2 随机游走。实际上没有必要这样做。我们只关心 t-score 而不是随机游走的那些实现。

对于这类问题，我们的目标是多次复制一个过程，首先编写一个函数来执行一次这样的过程是很方便的。你已经有了很好的工作代码；我们只需要将它包装在一个函数中（删除那些不必要的部分，如plot）：

sim <- function () {
  y <- cumsum(rnorm(100,0,1))
  x <- cumsum(rnorm(100,0,1))
  coef(summary(lm(y ~ x)))[2,3]
  }

这个函数不需要输入；它只返回一个实验的 t-score。

现在，我们将重复此操作 1000 次。我们可以编写一个for 循环，但函数replicate 更容易（如果需要，请阅读?replicate）

S <- replicate(1000, sim())

注意这需要一些时间，比这样一个简单的任务要慢得多，因为lm 和summary.lm 都很慢。稍后会展示一个更快的方法。

现在，S 是具有 1000 个值的向量，这就是您想要的“1000 个 t 检验的序列”。要获得“1,000 t-tests 的绝对值 > 1.96”的次数，我们可以使用

sum(abs(S) > 1.96)
# [1] 756

结果 756 正是我得到的；你会得到不同的东西，因为模拟是随机的。但正如预期的那样，它总是会很大。

更快的版本 sim：

fast_sim <- function () {
  y <- cumsum(rnorm(100,0,1))
  x <- cumsum(rnorm(100,0,1))
  y <- y - mean(y)
  x <- x - mean(x)
  xty <- crossprod(x,y)[1]
  xtx <- crossprod(x)[1]
  b <- xty / xtx
  sigma <- sqrt(sum((y - x * b) ^ 2) / 98)
  b * sqrt(xtx) * sigma
  }

这个函数计算没有lm的简单线性回归，以及没有summary.lm的t-score。

S <- replicate(1000, fast_sim())
sum(abs(S) > 1.96)
# [1] 778

另一种方法是使用cor.test：

fast_sim2 <- function () {
  y <- cumsum(rnorm(100,0,1))
  x <- cumsum(rnorm(100,0,1))
  unname(cor.test(x, y)[[1]])
  }

S <- replicate(1000, fast_sim())
sum(abs(S) > 1.96)
# [1] 775

让我们有一个基准：

system.time(replicate(1000, sim()))
#   user  system elapsed 
#  1.860   0.004   1.867 

system.time(replicate(1000, fast_sim()))
#   user  system elapsed 
#  0.088   0.000   0.090 

system.time(replicate(1000, fast_sim2()))
#   user  system elapsed 
#  0.308   0.004   0.312

cor.test 比lm + summary.lm 快得多，但手动计算更快！

【讨论】：