Monte Carlo模拟两个布朗运动的相关性(连续随机游走)
Monte Carlo simulation of correlation between two Brownian motion (continuous random walk)
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次(Monte Carlo次模拟),即围绕前面的程序建立一个循环,每次保存t统计。您将拥有一系列 1;000 个 t 检验:S = (t-test1, t-test2, ... , t-test1000)。计算 1,000 个 t 检验的绝对值 > 1.96 的次数,临界值在 5% 的显着性水平。如果该系列是 I(0),您会发现大约 5%。这里不会是这种情况(虚假回归)。
我需要添加什么来保存各自的系数?
您发布的与 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]
}
这个函数不需要输入;它只是 returns 一项实验的 t 分数。
现在,我们要重复这 1000 次。我们可以写一个 for
循环,但是函数 replicate
更容易(如有必要请阅读 ?replicate
)
S <- replicate(1000, sim())
请注意,这将需要一些时间,比完成如此简单的任务要慢得多,因为 lm
和 summary.lm
都很慢。稍后将展示一个更快的方法。
现在,S
是一个有1000个值的向量,就是你想要的"a sequence of 1000 t-tests"。要获得"the number of time the absolute value of the 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
快很多,但是人工计算更快!
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次(Monte Carlo次模拟),即围绕前面的程序建立一个循环,每次保存t统计。您将拥有一系列 1;000 个 t 检验:S = (t-test1, t-test2, ... , t-test1000)。计算 1,000 个 t 检验的绝对值 > 1.96 的次数,临界值在 5% 的显着性水平。如果该系列是 I(0),您会发现大约 5%。这里不会是这种情况(虚假回归)。
我需要添加什么来保存各自的系数?
您发布的与 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]
}
这个函数不需要输入;它只是 returns 一项实验的 t 分数。
现在,我们要重复这 1000 次。我们可以写一个 for
循环,但是函数 replicate
更容易(如有必要请阅读 ?replicate
)
S <- replicate(1000, sim())
请注意,这将需要一些时间,比完成如此简单的任务要慢得多,因为 lm
和 summary.lm
都很慢。稍后将展示一个更快的方法。
现在,S
是一个有1000个值的向量,就是你想要的"a sequence of 1000 t-tests"。要获得"the number of time the absolute value of the 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
.
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
快很多,但是人工计算更快!