Your script is rather inefficient with spurious cbind calls. Any particular reason not to use ?ar directly ?
Call: ar.yw.default(x = simtimeseries, order.max = 4) Coefficients: 1 2 3 4 1.9440 -1.9529 0.8450 -0.2154 Order selected 4 sigma^2 estimated as 15.29 To repeat the sim, you could use a for() loop but ?sapply is better: out<- sapply(1:100,function(...){ simtimeseries <- arima.sim(n=1024,list(order=c(4,0,0), ar=c(2.7607, -3.8106, 2.6535, -0.9258),sd=sqrt(1))) aryule <- ar.yw(simtimeseries,order.max=4) c( c(aryule$ar,NA)[1:4] , aryule$var.pred ) } ) rowMeans(out[1:4,]) # mean phi(1),...,4 see ?rowMeans for dealing with NA's mean(out[5,]) # mean sig^2 Cheers On Fri, Feb 10, 2012 at 6:42 AM, Jaymin Shah <jayminsh...@hotmail.com> wrote: > I have coded a time series from simulated data: > > simtimeseries <- arima.sim(n=1024,list(order=c(4,0,0),ar=c(2.7607, -3.8106, > 2.6535, -0.9258),sd=sqrt(1))) > #show roots are outside unit circle > plot.ts(simtimeseries, xlab="", ylab="", main="Time Series of Simulated Data") > > # Yule > ---------------------------------------------------------------------------- > > q1 <- cbind(simtimeseries[1:1024]) > q2 <- t(q1)%*%q1 > s0 <- q2/1204 > r1 <- cbind(simtimeseries[1:1023]) > r2 <- cbind(simtimeseries[2:1024]) > r3 <- t(r1)%*%r2 > s1 <- r3/1204 > t1 <- cbind(simtimeseries[1:1022]) > t2 <- cbind(simtimeseries[3:1024]) > t3 <- t(t1)%*%t2 > s2 <- t3/1204 > u1 <- cbind(simtimeseries[1:1021]) > u2 <- cbind(simtimeseries[4:1024]) > u3 <- t(u1)%*%u2 > s3 <- u3/1204 > v1 <- cbind(simtimeseries[1:1020]) > v2 <- cbind(simtimeseries[5:1024]) > v3 <- t(v1)%*%v2 > s4 <- v3/1204 > > i0 <- c(s0,s1,s2,s3) > i1 <- c(s1,s0,s1,s2) > i2 <- c(s2,s1,s0,s1) > i3 <- c(s3,s2,s1,s0) > > gamma <- cbind(i0,i1,i2,i3) > eta <-c(s1,s2,s3,s4) > inversegamma <- solve(gamma) > phihat <- inversegamma%*%eta > phihat > > Phihat <- cbind(phihat) > s <- c(s1,s2,s3,s4) > S <- cbind(s) > sigmasquaredyule <- s0 - (t(Phihat)%*%S) > sigmasquaredyule > > > > I did a yule walker estimate on the simulated data and wanted to work out phi > hat which is a vector of 4 values and sigmasquaredyule which is one value. > However, I want to run the simulated data 100 times i.e. in a for loop and > then take the averages of the phi hat values and sigmasquaredyule value. > > How would i repeat this simulated time series lots of times (e.g. a 100 > times) and store the average value of phi hat and sigmasquaredyule. > > Thank you > > > [[alternative HTML version deleted]] > > ______________________________________________ > 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. ______________________________________________ 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.