# ---------------------------------------------------------------------
# Program: LLA_Example1.S  
#  Author: Steven M. Boker
#    Date: June 17, 2002
#
# This program reads a simulated data file generated by the Mathematica
# program called DifferentialEQs1.nb.
# The simulated data is a single column in test.dat.
#
# Random measurement error is added to the simulated data (you'll want to
# remove that part).
#
# We then calculate derivatives with Local Linear Approximation and estimate
# the relationships between the derivatives
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- June 17, 2002
#      Create LLA_Example1.S
#
# ---------------------------------------------------------------------

options(width=100)

# ---------------------------------------------------------------------
# Read in the data from a single column textfile.

simData <- as.matrix(read.table("testLLA.dat"))

# ---------------------------------------------------------------------
# Here we add some simulated measurement error.  You'll want to remove this.

simError <- rnorm(length(simData),mean=0,sd=(0.5))

theData <- simData + simError

# ---------------------------------------------------------------------
# Now we set up the parameters for the analysis.

theInterval <- 1
minTau <- 1
maxTau <- 16

theTaus <- seq(minTau, maxTau, by=1)
totalN <- length(theData[,1])

# ---------------------------------------------------------------------
# Here we create a matrix to store the results of the analysis.

resultsMatrixLLA1 <- matrix(NA, length(theTaus), 9)

# ---------------------------------------------------------------------
# Calculate summary statistics and remove the linear trend and intercept.


tFrame <- data.frame(theData)
tLM <- lm(V1 ~ seq(1,(totalN * theInterval), length=totalN),data=tFrame)
summary(tLM)
theIntercept <- tLM$coef[1]
theTrend <- tLM$coef[2]
theResiduals <- tLM$residuals

# ---------------------------------------------------------------------
# Loop over the range of taus and fit the model.

j <- 1
for (theTau in theTaus) {
    t1 <- theResiduals[1:(totalN-(theTau*2))]
    t2 <- theResiduals[(1+theTau):(totalN-theTau)]
    t3 <- theResiduals[(1+(2* theTau)): totalN]
    tSelect <- !is.na(t1) & !is.na(t2) & !is.na(t3)
    tdataframe <- data.frame(
                    X=t2[tSelect], 
                    dX=(t3[tSelect]-t1[tSelect]) / (2*theTau*theInterval), 
                    d2X=(t3[tSelect] - (2 * t2[tSelect]) + t1[tSelect]) / (theTau^2 * theInterval^2))
    tLM <- lm(d2X ~ dX + X - 1, data=tdataframe)
    tVar <- var(tdataframe)
    resultsMatrixLLA1[j,1] <- theTau
    resultsMatrixLLA1[j,2] <- length(t2[tSelect])
    resultsMatrixLLA1[j,3] <- tLM$coefficients[1]
    resultsMatrixLLA1[j,4] <- tLM$coefficients[2]
    resultsMatrixLLA1[j,5] <- tVar[1,1]
    resultsMatrixLLA1[j,6] <- tVar[2,2]
    resultsMatrixLLA1[j,7] <- tVar[3,3]
    resultsMatrixLLA1[j,8] <- tVar[1,2]
    resultsMatrixLLA1[j,9] <- 1 - var(tLM$residuals) / var(tdataframe$d2X)
    j <- j + 1
}


colNames <- c("Tau", "N", "EstZeta", "EstEta", "VarX", "VarDX", "VarD2X", "CovXdX", "R^2")
rowNames <- rep(NULL, dim(resultsMatrixLLA1)[1])
dimnames(resultsMatrixLLA1) <- list(rowNames,colNames)

round(resultsMatrixLLA1,3)

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

# q()