############################################################################
# Generate some AR and MA sequences to see what they look like. Do realized
# spectra look like theoretical spectra? Can we tell models apart?
# Here are some MA(1) series: 4 series each for b = {-.8, -.2, .2, .8}
coeffs _ c ( -.8, -.2, .2, .8 )
nseries _ 4
n _ 100
MA _ array ( NA, c(length(coeffs),4,n) )
# Create the series
for ( i in seq(along=coeffs) )
for ( j in 1:nseries )
MA[i,j,] _ arima.sim ( model=list(ma=-coeffs[i]), n=100 )
# Plot the series
par ( mfrow=c(4,4), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) )
for ( j in 1:nseries )
ts.plot ( MA[i,j,], main=paste("beta =", coeffs[i]) )
mtext ( "Realizations of MA(1) series", outer=T, cex=2 )
# Plot the autocorrelations
par ( mfrow=c(4,5), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) ) {
# plot the theoretical autocorrelations
plot ( 0:15, c ( 1, coeffs[i]/(1+coeffs[i]^2), rep(0,14) ),
main="theoretical")
# plot the estimated autocorrelations
for ( j in 1:nseries )
acf ( MA[i,j,], lag.max=15 )
}
mtext ( "Autocorrelations of MA(1) series", outer=T, cex=2 )
# Plot the periodograms
par ( mfrow=c(4,5), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) ) {
# plot the spectrum
x _ seq ( 0, pi, length=40 )
plot ( x, 1 + 2*coeffs[i]*cos(x) + coeffs[i]^2, main="spectrum" )
# plot the periodogram
for ( j in 1:nseries )
spectrum ( MA[i,j,] )
}
mtext ( "Periodograms of MA(1) series", outer=T, cex=2 )
# Here are some AR(1) series: 4 series each for b = {-.8, -.2, .2, .8}
coeffs _ c ( -.8, -.2, .2, .8 )
nseries _ 4
n _ 100
AR _ array ( NA, c(length(coeffs),4,n) )
# Create the series
for ( i in seq(along=coeffs) )
for ( j in 1:nseries )
AR[i,j,] _ arima.sim ( model=list(ar=coeffs[i]), n=100 )
# Plot the series
par ( mfrow=c(4,4), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) )
for ( j in 1:nseries )
ts.plot ( AR[i,j,], main=paste("alpha =", coeffs[i]) )
mtext ( "Realizations of AR(1) series", outer=T, cex=2 )
# Plot the autocorrelations
par ( mfrow=c(4,5), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) ) {
# plot the theoretical autocorrelations
plot ( 0:15, c ( coeffs[i]^(0:15) ),
main="theoretical")
# plot the estimated autocorrelations
for ( j in 1:nseries )
acf ( AR[i,j,], lag.max=15 )
}
mtext ( "Autocorrelations of AR(1) series", outer=T, cex=2 )
# Plot the periodograms
par ( mfrow=c(4,5), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) ) {
# plot the spectrum
x _ seq ( 0, pi, length=40 )
plot ( x, ( (1-coeffs[i]*cos(x))^2 + (coeffs[i]*sin(x))^2 )^(-1),
main="spectrum" )
# plot the periodogram
for ( j in 1:nseries )
spectrum ( AR[i,j,] )
}
mtext ( "Periodograms of AR(1) series", outer=T, cex=2 )
# Plot the periodograms, based on an AR model
par ( mfrow=c(4,5), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=coeffs) ) {
# plot the spectrum
x _ seq ( 0, pi, length=40 )
plot ( x, ( (1-coeffs[i]*cos(x))^2 + (coeffs[i]*sin(x))^2 )^(-1),
main="spectrum" )
# plot the periodogram
for ( j in 1:nseries )
spectrum ( AR[i,j,], method="ar" )
}
mtext ( "Periodograms of AR(1) series", outer=T, cex=2 )
##########################################################################
# Now some AR(2) series to show some interesting behavior
n _ 100
a1 _ c ( .5, .2, -.2, -.5 )
a2 _ c ( -.2, -.5 )
AR2 _ array ( NA, c ( length(a1), length(a2), n ) )
# Create the series
for ( i in seq(along=a1) )
for ( j in seq(along=a2) ) {
print(c(i,j))
AR2[i,j,] _ arima.sim ( model=list(ar=c(a1[i],a2[j])), n=n )
}
# Plot the series
par ( mfrow=c(length(a1),length(a2)), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=a1) )
for ( j in seq(along=a2) )
ts.plot ( AR2[i,j,], main=paste("alpha =", a1[i], a2[j]) )
mtext ( "Realizations of AR(2) series", outer=T, cex=2 )
# Plot the autocorrelations
par ( mfrow=c(length(a1),length(a2)), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=a1) ) {
for ( j in seq(along=a2) ) {
# plot the estimated autocorrelations
acf ( AR[i,j,], lag.max=15 )
}
}
mtext ( "Autocorrelations of AR(1) series", outer=T, cex=2 )
# Plot the spectrum and periodogram
par ( mfrow=c(length(a1),length(a2)), oma=c(5.2,0,5.2,0) )
for ( i in seq(along=a1) ) {
for ( j in seq(along=a2) ) {
# plot the periodogram
spectrum ( AR2[i,j,], method="ar" )
# add the spectrum
x _ seq ( 0, pi, length=40 )
points ( x/(2*pi), ( (1-a1[i]*cos(x)-a2[j]*cos(2*x))^2 +
(a1[i]*sin(x)+a2[j]*sin(2*x))^2 )^(-1) )
}
}
mtext ( "Periodograms of AR(1) series", outer=T, cex=2 )