# Examples for the review of Fan and Yao # Synthetic data set # Logistic map corrupted with observational noise logistic.map <- function(x,r=4) { r*x*(1-x) } logistic.iteration <- function(n,x.init,r=4){ x <- vector(length=n) x[1] <- x.init for (i in 1:(n-1)) { x[i+1] <- logistic.map(x[i],r=r) } return(x) } x <- logistic.iteration(1000,x.init=runif(1)) y <- x+rnorm(1000,mean=0,sd=0.05) library(tseries) # Fit an autoregressive additive model # Uses spline smoothing, rather than local-linear regression with kernel # weights (which is what the review discusses and Fan and Yao favor), but # the differences are small, and I write less code this way # Inputs: time series (x), order of autoregression # Output: fitted GAM object am <- function(x,order) { require(mgcv) # Automatically generate a suitable data frame from the time series # and a formula to go along with it fit <- gam(as.formula(auto.formula(order)), data=design.matrix.from.ts(x,order)) return(fit) } # Convert a time series into a data frame of lagged values # Input: a time series, maximum lag to use # Output: a data frame with (order+1) columns, named lag0, lag1, ... , and # length(ts)-order rows design.matrix.from.ts <- function(ts,order) { n <- length(ts) x <- ts[(order+1):n] for (lag in 1:order) { x <- cbind(x,ts[(order+1-lag):(n-lag)]) } colnames(x) <- c("lag0",paste("lag",1:order,sep="")) return(as.data.frame(x)) } # Generate formula for an autoregressive GAM of a specified order # Input: order # Output: a formula which looks like # "lag0 ~ s(lag1) + s(lag2) + ... + s(lagorder)" auto.formula <- function(order) { inputs <- paste("s(lag",1:order,")",sep="",collapse="+") form <- paste("lag0 ~ ",inputs) return(form) } # Plot the first part of the time series plot(y[1:100],xlab="t",ylab=expression(y[t]),type="l") # Plot successive values of y against each other plot(lag0 ~ lag1,data=design.matrix.from.ts(y,1),xlab=expression(y[t]), ylab=expression(y[t+1]),pch=16) # Add the linear regression (which would be the AR(1) model) abline(lm(lag0~lag1,data=design.matrix.from.ts(y,1)),col="red") # Fit an AR(8) and add its fitted values yar8 <- arma(y,order=c(8,0)) points(y,fitted(yar8),col="red") # Fit a first-order nonparametric autoregression, add fitted values yaar1 <- am(y,order=1) points(y[-length(y)],fitted(yaar1),col="blue")