我想以下的回归做10000模拟,以便找到H:模拟R中回归
日志(TV)= H *日志(M)
在日志(TV)是布朗运动总体变化的对数(BM),m是x轴上不相交子区间的大小。所以,通过改变m的大小,我们对电视和日志(电视)具有新的价值。运行日志(TV)对log(m),我们可以找到H作为拟合线的斜率。以下是我一直在努力代码:
t<-0:1000
nsim <- 10000
sig2<-0.01
#simulate with single BM
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
d1<- c(0, cumsum(x))
# calculating total variation
v1<-d1[seq(1,length(d1),1)] # m=1
v2<-d1[seq(1,length(d1),10)] # m=10
v3<-d1[seq(1,length(d1),100)]#m=100
s1<-sum(abs(v1[1:(length(v1)-1)]-v1[2:length(v1)]))
s2<-sum(abs(v2[1:(length(v2)-1)]-v2[2:length(v2)]))
s3<-sum(abs(v3[1:(length(v3)-1)]-v3[2:length(v3)]))
s<-c(log(s1),log(s2),log(s3))
m<-c(log(1),log(10),log(100))
fit<-lm(s~m)
这可以帮助我找到h的一个值,现在我想做10000个模拟和下面是我想出了代码:
for (i in 1:nsim) {X[i, ] <- c(0, cumsum(rnorm(n = length(t) - 1, sd = sqrt(sig2))))}#10000 simulations of BM
# calculating total variation:
v1<-X[,1:1==1]
v2<-X[,1:3==10]
v3<-X[,1:100==100]
a1<-abs(v1[,1:ncol(v1)-1]-v1[,2:ncol(v1)])
a2<-abs(v2[,1:ncol(v2)-1]-v2[,2:ncol(v2)])
a3<-abs(v3[,1:ncol(v3)-1]-v1[,2:ncol(v3)])
s1<-apply(a1,1,sum)
s2<-apply(a2,1,sum)
s3<-apply(a3,1,sum)
s<-cbind(s1,s2,s3)
S<-log(s)
M<-c(log(1),log(10),log(100))
矩阵S是其中每行将包含3个对应于m = 1,10,100的log(TV)值的矩阵。现在我需要在S中用M对每一行进行回归,并在每次运行回归时记录H的值。我被困在这一步。有人可以建议我的方式,我可以做到这一点?
非常感谢你提前任何建议
的'replicate'功能是一个标准的方法这个。 's'和'm'只有长度三的事实不是你的担心吗? – 2014-11-21 00:25:57