2014-11-20 90 views
0

我想以下的迴歸做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的值。我被困在這一步。有人可以建議我的方式,我可以做到這一點?

非常感謝你提前任何建議

+1

的'replicate'功能是一個標準的方法這個。 's'和'm'只有長度三的事實不是你的擔心嗎? – 2014-11-21 00:25:57

回答

2

,說明如何與原來的代碼中使用複製而不是什麼在你的循環:

Res <- replicate(n= 100, expr = { 
t<-0:1000 
nsim <- 10000 
sig2<-0.01 

x <- rnorm(n = length(t) - 1, sd = sqrt(sig2)) 
d1<- c(0, cumsum(x)) 

v1<-d1[seq(1,length(d1),1)] 
v2<-d1[seq(1,length(d1),10)] 
v3<-d1[seq(1,length(d1),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); coef(fit) 
    } 
) 

> str(Res) 
num [1:2, 1:100] 4.299 -0.526 4.328 -0.506 4.441 ... 
- attr(*, "dimnames")=List of 2 
    ..$ : chr [1:2] "(Intercept)" "m" 
    ..$ : NULL 
+0

非常感謝您的建議。這是完美的,關於你對s和m的長度的關注,我試圖找到一種方法,我可以概括生成v的過程。我使用循環,但它似乎不能很好地工作 – Ocean 2014-11-21 01:56:18