[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