[R] help on R for loops

From: song xin <justinxin_at_yahoo.com>
Date: Sat, 26 Jan 2008 21:24:20 -0800 (PST)


Hi, all.
  I need help on improving the efficiency of my R simulations. Below is a function that simulates a change point model. It first generates a sequence of three dimensional ARMA(1,1) observations, then calculates the one step ahead prediction errors, some statistic is calcualted and compared with threshold values in the end. As you can see in the function, there are 5 for loops, which makes the simulation long and inefficient. Can anyone help me improve it?

simfun<-function(b){
for(l in 1:20)
{
qqqqqq<-matrix(0,2,40)
for (w in 1:40)
{
x2 <- matrix(rep(0,6000),nr=3)
epsilonzero2<-(mvrnorm(n=1,mu=c(0,0,0),Sigma=residsigma)) x2[,1:2]<-c(0,0,0)
for (i in 3:550)
{
epsilonone<-(mvrnorm(n=1,mu=c(0,0,0),Sigma=residsigma))   x2[,i] <-
phi%*%x2[,i-1]+epsilonone-theta%*%epsilonzero2 epsilonzero2<-epsilonone
}

residsigma3<-matrix(c(0.789,0.2143,0.171,0.2143,1.4394,-0.229,0.171,-0.229,0.6649),nrow=3,byrow=T) epsilonzero3<-epsilonzero2
for (i in 551:2000)
{
epsilonone<-(mvrnorm(n=1,mu=c(0,0,0),Sigma=residsigma3))   x2[,i] <-
phi%*%x2[,(i-1)]+epsilonone-theta%*%epsilonzero3 epsilonzero3<-epsilonone
}

inno2<-matrix(0,3,2000)
inno2[,1]<-c(0,0,1)
for ( i in 2:2000)
{
inno2[,i]<-x2[,i]-phi%*%x2[,(i-1)]+
theta%*%inno2[,(i-1)]
}

sampeigen4<-matrix(0,3,(2000-nnumber[l])) for ( i in (nnumber[l]+1):2000)
{
var1<-matrix(0,3,3)
for ( j in 1:(nnumber[l]))
{
var1<-var1+j^{b}*inno2[,i-nnumber[l]-1+j]%*%t(inno2[,i-nnumber[l]-1+j])/(sum(c(1:nnumber[l])^{b})) var1<-var1
}

sampeigen4[,(i-nnumber[l])]<-(eigen(var1)$values-eigen(residsigma)$values)
}

chisqstat<-rep(0,(2000-nnumber[l]))
for(i in 1:(2000-nnumber[l]))
{
chisqstat[i]<-((nnumber[l]-1)/2)*t(sampeigen4[,i])%*%(diag(eigen(residsigma)$values)^2)%*%(sampeigen4[,i])
}

normstat<-apply((sampeigen4+eigen(residsigma)$values),2,prod) qqqqqq[1,w]<-vp(chisqstat[(550-nnumber[l]):(2000-nnumber[l])]) qqqqqq[2,w]<-vp1(normstat[(550-nnumber[l]):(2000-nnumber[l])])
}

inarl2[,l]<-apply(qqqqqq,1,mean)
}

inarl2
}

Thanks

XS



Looking for last minute shopping deals?

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 27 Jan 2008 - 05:28:49 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 27 Jan 2008 - 05:30:11 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive