# Simulating Martindale's so-called analysis of
# "primordial content" in Beethoven's music
  # After pp. 286--289 of Martindale's _The Clockwork Muse_


require(zoo) # For moving average function


# Generate white noise and subject it to Martindale's procedures
# Inputs: flag for whether to make a plot or not
# Outputs: returns the fitted quadratic-and-AR(2) model
rmartindale <- function(make.plot=TRUE) {
  # Generate a completely random series
    # One year longer than what Martindale displays, because he gives a 2-year
    # moving average
  years <- 1794:1826
  n <- length(years)
  x <- rnorm(n)
  # Take the moving average with a 2-year window
  y = rollmean(x,2)
  # Fit a quadratic on year to the moving average
  quad <- lm(y ~ poly(years[-1],2))
    # Use orthogonal polynomials for numerical stability

  # Following definitions not really necessary but ease comparison with
  # Martindale's results:
  t <- years-1794
  t.sq <- t^2
  # Fit a quadratic time trend plus second-order autoregression
  quad.and.ar <- lm(y[-(1:2)] ~ t[-(1:3)]+ t.sq[-(1:3)] + y[-c(1,n-1)] + y[-((n-2):(n-1))])

  if (make.plot) {
    plot(years[-1],y,type="b", xlab="year", ylab="(pseudo-) Primordial Content")
    lines(years[-1],fitted(quad),lty=2)
  }
  return(quad.and.ar)
}

# Get the sampling distribution of R^2 when Martindale's procedure is applied
# to Gaussian white noise
  # Given as a function because running it for large numbers of replicates can
  # be time-consuming; this way you can source the file without having it
  # lock up R
# Inputs: Number of replicates B
# Output: Vector of R^2 values of length B
# Calls: rmartindale, of course
sampling.dist.of.r.squared <- function(B) {
  return(replicate(B,summary(rmartindale(make.plot=FALSE))$r.squared))
}

# Get the sampling distribution of AR(2) coefficients from moving averages of
# white noise
  # Including the effects of fitting a quadratic trend simultaneously
# Inputs: Number of replciates B
# Outputs: 2*B matrix of AR coefficients
sampling.dist.of.ar.coefs <- function(B) {
  return(replicate(B,summary(rmartindale(make.plot=FALSE))$coefficients[4:5]))
}