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@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.
