## This function produces the innovariogram/regressogram, as in Figure 4.8.
## It calls another function, Stats.ar
innovariogram.regressogram<-function(data, MSE){
  V <- Stats.ar(data, MSE) 
  n <- (ncol(data)-1)/2       #for time index
  m <- ncol(data)-1           #for lag index
  par(mfrow=c(2:1))      
  par(pch=16)

  #sample log innovation variance
  logdelta <- log(diag(V)) 
  # prepare an empty plot with correct limits and no axes

  plot(range(c(-n,n)),range(logdelta[!is.na(logdelta)])+c(-0.7,0.4),
  type = "n", axes = T, xlab = "Time index", ylab = "Log innovation variance")  

  # show points for log innovation variance vs time index
  points(as.vector(c(-n:n)),as.vector(logdelta),cex=0.5)

  # sample autoregressive coefficient
  autoreg <- sapply(2:(m+1), function(i) if(i==(m+1)) V[i,1] else diag(V[i:(m+1),1:(m+2-i)]))
    
  # prepare an empty plot with correct limits and no axes
  plot(range(c(0,m)),range(autoreg[!is.na(autoreg)])+c(-0.2,0.1),
  type = "n", axes = T, xlab = "Lag", ylab = "Autoregressive coefficient")

  # show points for log autoregressive coefficient vs lag
  for(i in 2:ncol(data)){
    coef <- autoreg[[i-1]]
    points(as.vector(c(rep(i-1, ncol(data)-i+1))), as.vector(coef), cex=0.5)
  }
}