# Produces chart like Chart C in G E P Box and G M Jenkins, Time Series
# Analysis: Forecasting and Control, San Francisco: Holden-Day 1970 and 1976.
# 
par(col=1)
plot(0,type="n",xlim=c(-2,2),xlab="theta1",ylim=c(-1,1),ylab="theta2",
     axes=T)
for (i in 1:2) segments(-6+4*i,-1,0,1)
k <- -0.47
par(col=2,new=T)
curve(0.5*(-1/k-sqrt(1/k^2-4*(1+x^2))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)
for (i in 1:4){
  k <- -i/10
  par(new=T)
  curve(0.5*(-1/k-sqrt(1/k^2-4*(1+x^2))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
}
abline(h=0)
for (i in 1:4){
  k <- i/10
  par(new=T)
  curve(0.5*(-1/k+sqrt(1/k^2-4*(1+x^2))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
}
k <- 0.47
par(new=T)
curve(0.5*(-1/k+sqrt(1/k^2-4*(1+x^2))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)
par(col=3)
for (i in 1:6){
  par(new=T)
  k <- -i/10
  curve(0.5*(x/k+sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
  par(new=T)
  curve(0.5*(x/k-sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
}
par(new=T)
k <- -2/3
curve(0.5*(x/k+sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)
par(new=T)
curve(0.5*(x/k-sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)
abline(v=0)
for (i in 1:6){
  par(new=T)
  k <- i/10
  curve(0.5*(x/k+sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
  par(new=T)
  curve(0.5*(x/k-sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
    ylim=c(-1,1),ylab="",axes=F)
}
par(new=T)
k <- 2/3
curve(0.5*(x/k+sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)
par(new=T)
curve(0.5*(x/k-sqrt(x^2/k^2-4*(1+x^2+x/k))),-2,2,xlim=c(-2,2),xlab="",
  ylim=c(-1,1),ylab="",axes=F)