[R] help on R for loops

song xin justinxin at yahoo.com
Sun Jan 27 06:24:20 CET 2008


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?



More information about the R-help mailing list