[R] Help needed please

[ilai]

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 <jayminshah1 at 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 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.