# ------------------------------------------------------
# Program: TimeSeries3.S  
#  Author: Steven M. Boker
#    Date: Tue Nov 23 10:18:26 EST 2004
#
#  Timeseries, State Space and Vector Plots Part 3
#
# ------------------------------------------------------

# ------------------------------------------------------
#  Simulate a variable that uses "Greater than 10" as a
#  category to chop off a long tail.

tX <- (rchisq(300, df=5)^3)/20
tY <- .8*tX + rnorm(300, mean=0, sd=2)
tXCeiling <- tX
tXCeiling[tX > 10] <- 10
tXCeiling <- round(tXCeiling,0)
tYCeiling <- tY
tYCeiling[tY > 2] <- 2

tFrame <- data.frame(tX, tY, tXCeiling)


# ------------------------------------------------------
#  Calculate a smoothed bivariate density function using a
#  Rectangular kernel for Simulated "Two Population Data"

smoothX <- seq(0, 10, by=1)
smoothY <- seq(-2, 17, by=2)
smoothF <- matrix(NA, length(smoothX), length(smoothY))
h <- 1
for (i in 1:length(smoothX)) {
   x <- smoothX[i]
   t1x <- abs(x - tXCeiling)/h
   for (j in 1:length(smoothY)) {
      y <- smoothY[j]
      t1y <- abs(y - tY)/h
      t2 <- rep(0, length(t1x))
      t2[t1x < 1 & t1y < 1 & !is.na(t1x) & !is.na(t1y)] <- 1/2
      smoothF[i,j] <- sum(t2)/(length(tXCeiling)*h) 
   }
}
smoothF[is.na(smoothF)] <- 0
smoothF <- smoothF/sum(smoothF)

# ------------------------------------------------------
#  Plot a line contour graph

graphsheet(height=6.4,width=6.4)
contour(smoothX, smoothY, smoothF, nlevels=8,
        xlab = "tXCeiling", 
        ylab = "tY")


tFrame <- data.frame(tXCeiling=rep(smoothX,length(smoothY)),
                     tY=rep(smoothY,each=length(smoothX)),
                     Prob=c(smoothF))



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

tX <- sin(seq(0,20, by=.1))
tY <- 5 + .5 * tX[1:201] + rnorm(201,mean=0,sd=.05)
tX <- tX[1:201] + rnorm(201,mean=0,sd=.05)

tData <- data.frame(x=tX,y=tY)
summary(tData)

# ------------------------------------------------------
#  Create a timelagged dataframe (3-D Embedding)

tau <- 15
tLen <- length(tX)
rDataLagged <- data.frame(x1=tX[1:(tLen-(2*tau))], 
                          x2=tX[(1+tau):(tLen-tau)],
                          x3=tX[(1+(2*tau)):tLen],
                          y1=tY[1:(tLen-(2*tau))], 
                          y2=tY[(1+tau):(tLen-tau)],
                          y3=tY[(1+(2*tau)):tLen])

# ------------------------------------------------------
#  Create an cross-correlation function plot of tX and tY
 
maxTau <- 50
tLen <- length(tX)
tCCF1 <- rep(NA, maxTau+1)

for (tau in (-maxTau:maxTau)) {
   t1 <- tX[(1+maxTau):(tLen-maxTau-tau)]
   t2 <- tY[(1+maxTau+tau):(tLen-maxTau)]
   tSelect <- !is.na(t1) & !is.na(t2)
   if (length(t1[tSelect]) < 5) next
   tCCF1[maxTau+tau+1] <- cor(t1[tSelect], t2[tSelect])
}

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


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

tX <- sin(seq(0,20, by=.1))
tY <- 5 + .5 * tX[21:201] + rnorm(181,mean=0,sd=.05)
tX <- tX[1:181] + rnorm(181,mean=0,sd=.05)

tData <- data.frame(x=tX,y=tY)
summary(tData)


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

write(tY, "c:/Class/Psych344509/CrossSine1.dat", 
      ncolumns=1)


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