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:
1234
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.