CODEWITHSUNDEEP
r
Ocean
Ocean

Reputation: 81

Simulation for regression in R

I am trying to do 10000 simulation for the following regression in order to find H:

log(TV)=H*log(m)

in which log(TV) is the log of total variation of Brownian motion (BM) and m is the size of the disjoint subinterval on x-axis. So, by changing size of m, we have new value for TV and thus log(TV). Running log(TV) against log(m) we can find H as the slope of the fitted line. The following is the code that I have been working on:

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)

This helps me to find one value of H and now I want to do 10000 simulations and the following is the code that I came up with:

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))

Matrix S is the matrix in which each row will contain 3 values of log(TV) that correspond to m=1,10,100. Now I need to run regression of each row in S with M and record the value of H each time I run regression. I got stuck at this step. Can someone please suggest me the way that I can do that??

Thank you so much in advance for any suggestions

Upvotes: 0

Views: 223

Answers (1)

IRTFM
IRTFM

Reputation: 263481

Illustrating how to use replicate with the original code rather than what was in your loop:

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

Upvotes: 2

Related Questions