# Plots autoregression coefficient estimates, the autocorrelation # function, and the partial autocorrelation function for an # autoregressive model on the residuals of a temperature model fit in # other temp.R files. # # Figure caption: Autoregressive model of order p = 22 for core # temperature residuals. TOP: Coefficients theta-hat-i as a function # of lag i. MIDDLE: the sample autocorrelation function. BOTTOM: # The sample partial autocorrelation function. # The Haell Miscellaneous library contains a variety of utility # functions for statistical analysis, data manipulation, and LaTeX # expressions. library(Hmisc) # First read in the temperature data, and assign column names. temperatureDF <- read.table("data/temp.dat", col.names = c("Temp", "Time")) temperature <- temperatureDF$Temp time <- temperatureDF$Time # We now want to regress our data against sinusoids with a period of 72. # First create the sinusoids. cosine <- cos(2 * pi * time/72) sine <- sin(2 * pi * time/72) # Then fit a linear model. lm.temperature <- lm(temperature ~ cosine + sine) temperatureResiduals <- lm.temperature$residuals orderP <- 22 n <- 330 x <- array(0, dim = c(n, orderP)) range <- (orderP + 1):(orderP + n) y <- temperatureResiduals[range] for(i in 1:22) { x[ , i] <- Lag(temperatureResiduals, i)[range] } # Fit autoregressive model without intercept. ar.reg=lm(y ~ x - 1) summary(ar.reg) # Display summary. # Set graphical parameters. par(oma = rep(2, 4)) par(mar = c(5, 5, 1, 1)) par(mfrow=c(3,1)) axis.size <- 1.4 label.size <- 1.8 # Plot estimated autocorrelation coefficients. plot(1:22, ar.reg$coef, type="l", xlab = "Lag", ylab = expression(phi[i]), lwd=2, cex.lab = label.size, cex.axis = axis.size, yaxt = "n") abline(h = 0, lty = 3, lwd = 2) axis(2, at = seq(from = -0.2, to = 1, by = 0.4), labels = seq(from = -0.2, to = 1, by = 0.4), cex.axis = 1.5) # Plot sample autocorrelation function. acf(temperatureResiduals, main = "", cex.lab = label.size, cex.axis = axis.size, lwd = 2, yaxt = "n", ci.col = 1) axis(2, at = seq(from = 0, to = 1, by = 0.5), labels = seq(from = 0, to = 1, by = 0.5), cex.axis = axis.size) # Plot sample partial autocorrelation function. pacf(temperatureResiduals, main = "", cex.lab = label.size, cex.axis = axis.size, lwd = 2, yaxt = "n", ci.col = 1) axis(2, at = seq(from = -0.3, to = 0.6, by = 0.3), labels = seq(from = -0.3, to = 0.6, by = 0.3), cex.axis = axis.size) dev.print(device = postscript, "18.7.eps", horizontal = TRUE)