# ------------------------------------------------------
# Program: TimeSeries1.S  
#  Author: Steven M. Boker
#    Date: Tue Nov 16 12:37:33 EST 2004
#
#  Timeseries Plots Part 1
#
# ------------------------------------------------------


# ------------------------------------------------------
#  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 <- 5 + .5 * tX1 + rnorm(201, mean=0, sd=5)
tX1 <- tX1 + rnorm(201, mean=0, sd=5)
tX1[tX1 < 0] <- 0

tData1 <- data.frame(x=tX1, y=tY1)
summary(tData1)

# ------------------------------------------------------
#  Explore Transforms of tX1
 
t1 <- tX1

t1 <- log(tX1)

# ------------------------------------------------------
#  Create a timeseries of tX1
 
tTime <- seq(0,1000, by=5)/60
tData1a <- data.frame(x=tX1, y=tY1, time=tTime)

# ------------------------------------------------------
#  Create a scrambled timeseries of tX1
 
tTimeRand <- order(runif(201,0,1))
tData1b <- data.frame(x=tX1, y=tY1, time=tTimeRand)

# ------------------------------------------------------
#  Create a time-lagged timeseries of tX1
 
tau <- 5
tLen <- length(tX1)
tDataLag1 <- data.frame(x1=tX1[1:(tLen-tau)], x2=tX1[(1+tau):tLen])

# ------------------------------------------------------
#  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=3)
}

tY2 <- 5 + .5 * tX2 + rnorm(201, mean=0, sd=5)
tX2 <- tX2 + rnorm(201, mean=0, sd=2)

tData2 <- data.frame(x=tX2, y=tY2)
summary(tData2)

# ------------------------------------------------------
#  Create a scrambled timeseries of tX2
 
tTimeRand2 <- order(runif(201,0,1))
tData2b <- data.frame(x=tX2[tTimeRand], time=1:201)

# ------------------------------------------------------
#  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])


# ------------------------------------------------------
#  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])
}

graphsheet(height=6.4,width=7.5)
plot(c(0,maxTau), c(-1,1), type="n" ,
   xlab="Lag",
   ylab="Autocorrelation")
lines(c(0:maxTau), tACF1, type="l", lty=1)
lines(c(0,maxTau), c(0,0), type="l", lty=2)


# ------------------------------------------------------
#  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])
}

graphsheet(height=6.4,width=7.5)
plot(c(0,maxTau), c(-1,1), type="n" ,
   xlab="Lag",
   ylab="Autocorrelation")
lines(c(0:maxTau), tACF2, type="l", lty=1)
lines(c(0,maxTau), c(0,0), type="l", lty=2)


# ------------------------------------------------------
#  Create a cyclic time series and 
#  synchronous cross-regression relation

tX3 <- sin(seq(0,20, by=.1))
tY3 <- 5 + .5 * tX3 + rnorm(201,mean=0,sd=1)
tX3 <- tX3 + rnorm(201,mean=0,sd=1)

tData3 <- data.frame(x=tX3,y=tY3)
summary(tData3)


# ------------------------------------------------------
#  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])
}

graphsheet(height=6.4,width=7.5)
plot(c(0,maxTau), c(-1,1), type="n" ,
   xlab="Lag",
   ylab="Autocorrelation")
lines(c(0:maxTau), tACF3, type="l", lty=1)
lines(c(0,maxTau), c(0,0), type="l", lty=2)

# ------------------------------------------------------
#  Create a time-lagged timeseries of tX2
 
tau <- floor(60/4)
tLen <- length(tX3)
tDataLag3 <- data.frame(x1=tX3[1:(tLen-tau)], x2=tX3[(1+tau):tLen])


# ------------------------------------------------------
#  Write the timeseries matrix to a textfile for input to
#  a recurrence analysis program.

tMatrix <- cbind(tX3, tY3)
write(round(t(tMatrix),5), 
      "c:/Class/Psych344509/CrossSine1.dat", 
      ncolumns=dim(tMatrix)[2])

# ------------------------------------------------------
# Quit the program