# ------------------------------------------------------ # Program: TimeSeries2.R # Author: Steven M. Boker # Date: Tue Oct 31 09:26:00 EST 2006 # # Timeseries, State Space and Vector Plots Part 2 # # ------------------------------------------------------ # ------------------------------------------------ # Load required libraries. library(lattice) # ------------------------------------------------------ # Create an error-free autoregressive time series and # synchronous cross-regression relation tX1 <- rep(NA, 201) tX1[1] <- 90 for (i in 2:201) { tX1[i] <- .975 * tX1[i-1] } tY1 <- 10 + .5 * tX1 + rnorm(201, mean=0, sd=5) tX1 <- 10 + tX1 + rnorm(201, mean=0, sd=5) tX1[tX1 <= 0] <- 0.001 # ------------------------------------------------------ # Create a timeseries of tX1 tTime <- seq(0,1000, by=5)/60 tData1a <- data.frame(x=tX1, y=tY1, time=tTime) # ------------------------------------------------------ # Create an error dependent autoregressive time series and # synchronous cross-regression relation tX2 <- rep(NA, 201) tX2[1] <- 1 for (i in 2:201) { tX2[i] <- .9 * tX2[i-1] + rnorm(1, mean=0, sd=2) } tY2 <- 5 + .5 * tX2 + rnorm(201, mean=0, sd=5) tX2 <- tX2 + rnorm(201, mean=0, sd=1) tTime <- seq(0,1000, by=5)/60 tData2 <- data.frame(x=tX2, y=tY2, time=tTime) summary(tData2) pdf("TimeSeries2TSPlot1.pdf", height=5, width=6) plot(tData1a$time, tData1a$x, xlab="Time", ylab="X", type='l', pch=1) dev.off() pdf("TimeSeries2TSPlot2.pdf", height=5, width=6) plot(tData2$time, tData2$x, xlab="Time", ylab="X", type='l', pch=1) dev.off() # ------------------------------------------------------ # Create a scrambled timeseries of tX2 tTimeRand2 <- order(runif(201,0,1)) tData2b <- data.frame(x=tX2[tTimeRand2], time=1:201) pdf("TimeSeries2TSPlot3.pdf", height=5, width=6) plot(tData2b$time, tData2b$x, xlab="Time", ylab="X", type='l', pch=1) dev.off() # ------------------------------------------------------ # Create a time-lagged timeseries of tX2 tau <- 3 tLen <- length(tX2) tDataLag2 <- data.frame(x1=tX2[1:(tLen-tau)], x2=tX2[(1+tau):tLen]) pdf("TimeSeries2SSPlot1.pdf", height=5, width=5) plot(tDataLag2$x1, tDataLag2$x2, xlab="X(t)", ylab="X(t+3)", main="Random Walk State Space", type='l', pch=1) dev.off() tau <- 3 tLen <- length(tX2) tDataLag2Scram <- data.frame(x1=tData2b$x[1:(tLen-tau)], x2=tData2b$x[(1+tau):tLen]) pdf("TimeSeries2SSPlot2.pdf", height=5, width=5) plot(tDataLag2Scram$x1, tDataLag2Scram$x2, xlab="X(t)", ylab="X(t+3)", main="Scrambled Random Walk State Space", type='l', pch=1) dev.off() # ------------------------------------------------------ # Create an autocorrelation function plot of tX1 maxTau <- 100 tLen <- length(tX1) tACF1 <- rep(NA, maxTau+1) for (tau in (0:maxTau)) { t1 <- tX1[1:(tLen-tau)] t2 <- tX1[(1+tau):tLen] tSelect <- !is.na(t1) & !is.na(t2) if (length(t1[tSelect]) < 5) next tACF1[tau+1] <- cor(t1[tSelect], t2[tSelect]) } backward <- seq(maxTau+1, 2, by= -1) pdf("TimeSeries2ACF1.pdf", height=5, width=6) plot(c(-maxTau,maxTau), c(-1,1), type="n", xlab="Lag", ylab="Autocorrelation") lines(c(-maxTau:maxTau), c(tACF1[backward],tACF1), type="l", lty=1) lines(c(-maxTau,maxTau), c(0,0), type="l", lty=2) dev.off() # ------------------------------------------------------ # Create an autocorrelation function plot of tX2 maxTau <- 100 tLen <- length(tX2) tACF2 <- rep(NA, maxTau+1) for (tau in (0:maxTau)) { t1 <- tX2[1:(tLen-tau)] t2 <- tX2[(1+tau):tLen] tSelect <- !is.na(t1) & !is.na(t2) if (length(t1[tSelect]) < 5) next tACF2[tau+1] <- cor(t1[tSelect], t2[tSelect]) } backward <- seq(maxTau+1, 2, by= -1) pdf("TimeSeries2ACF2.pdf", height=5, width=6) plot(c(-maxTau,maxTau), c(-1,1), type="n", xlab="Lag", ylab="Autocorrelation") lines(c(-maxTau:maxTau), c(tACF2[backward],tACF2), type="l", lty=1) lines(c(-maxTau,maxTau), c(0,0), type="l", lty=2) dev.off() # ------------------------------------------------------ # Create a cyclic time series and # synchronous cross-regression relation tX3 <- sin(seq(0,20, by=.2)) * .5 tY3 <- 5 + .5 * tX3 + rnorm(101,mean=0, sd=.5) tX3 <- tX3 + rnorm(101,mean=0,sd=.5) tData3 <- data.frame(x=tX3,y=tY3) summary(tData3) # ------------------------------------------------------ # Create a timeseries plot of tX3 tLen <- length(tX3) pdf("TimeSeries2TSPlot4.pdf", height=5, width=6) plot(c(0,tLen), c(-2,2), xlab="Time", ylab="X", type='n', pch=1) lines(c(1:tLen), tX3, type="l", lty=1) lines(c(0,tLen), c(0,0), type="l", lty=2) dev.off() # ------------------------------------------------------ # Create an autocorrelation function plot of tX3 maxTau <- 100 tLen <- length(tX3) tACF3 <- rep(NA, maxTau+1) for (tau in (0:maxTau)) { t1 <- tX3[1:(tLen-tau)] t2 <- tX3[(1+tau):tLen] tSelect <- !is.na(t1) & !is.na(t2) if (length(t1[tSelect]) < 5) next tACF3[tau+1] <- cor(t1[tSelect], t2[tSelect]) } backward <- seq(maxTau+1, 2, by= -1) pdf("TimeSeries2ACF3.pdf", height=5, width=6) plot(c(-maxTau,maxTau), c(-1,1), type="n", xlab="Lag", ylab="Autocorrelation") lines(c(-maxTau:maxTau), c(tACF3[backward],tACF3), type="l", lty=1) lines(c(-maxTau,maxTau), c(0,0), type="l", lty=2) dev.off() # ------------------------------------------------------ # Create a State Space Plot of tX3 tLen <- length(tX3) for(tTau in seq(1, 30, by=2)) { tFrame <- data.frame(x1=tX3[1:(tLen-tTau)], x2=tX3[(1+tTau):tLen]) tFile <- paste("TimeSeries2SSTau", tTau, ".pdf", sep="") pdf(tFile, height=5, width=5) plot(c(-3, 3), c(-3, 3), xlab="X(t)", ylab=paste("X(t+", tTau, ")", sep=""), main="Cyclic State Space", type='n') lines(tFrame$x1, tFrame$x2, pch=1, type='p') abline(lm(x2 ~ x1, data = tFrame)) dev.off() } # ------------------------------------------------------ # Write the timeseries matrix to a textfile for input to # a recurrence analysis program. tMatrix <- as.matrix(tData3) write(round(t(tMatrix),4)[1,], "tSineData1.dat", ncolumns=1) # ------------------------------------------------------ # Quit the program