#========================================================== # Simulation of an AR(1) Model using arima.sim() # Code is from Z&W (slightly modified) # Also illustrates par(mfrow(2,2)) and graph labeling #======================================================= #The simulation using start.innov burn-in. We simply add "something" (here it is 1) #to get a series that has a mean differing from zero. e<-rnorm(100,sd=1) e.start=rnorm(25, sd=1) y.ar1<- 1+arima.sim(model=list(ar=0.75), n=100,innov=e, start.innov=e.start) #Here the values in start.innov are used for the "burn in." For a nice AR(1) model #a burn-in period of 25 probably suffices, but for more complicated models you may want a #longer burn-in period. #By the way, you "need" just as much burn-in to generate a series of length 100 as to #generate one of length 5000. # Let's check on the mean to see if it makes sense. # We may run our simulation several times. mean(y.ar1) # If we want to have the simulation begin with the same random seed each time # we can explicitly set the pseudo-random number generator seed set.seed(101) # When we use this seed we get a larger mean than we might have expected, but this # but due to the reasonably large variance of this series (calculate tau!), the mean # will bounce a lot from seed to seed. #Now lets do some pictures --- and put three on just one page: # #(1) the simulated series, #(2) the theoretical ACF #(3) the empirical ACF for our simulation par(mfrow=c(2,2)) tsplot(y.ar1, main="Simulated AR(1)") abline(h=1) gamma.j<-rep(.75,10)^(0:9) tsplot(gamma.j, type="h", main="Theoretical ACF for AR(1)", ylab="Autocorreleation", xlab="Lag+1") temp<-acf(y.ar1, lag.max=10, plot=F) tsplot(temp\$acf,type="h", main="Empirical ACF for Simulated AR(1)", ylab="Autocorreleation", xlab="Lag+1") #One line to make it prettier: abline(h=0) #How do the series compare? To me the typical run looks pretty reasonable, #though you may notice that at higher lags the sample ACF tends to "turn up." #Thus, there is a sampling effect that is not present in the theoretical ACF #You might want to look at the help screens for the acf function.