# ---------------------------------------------------------------------
# Program: LDE_Example2.S  
#  Author: Steven M. Boker
#    Date: June 18, 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).
#
# The Mx script to run the model is created by the S code.  This way
# we can simulate across a variety of taus, number of indicators and methods
# for estimating the second derivative.
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- June 18, 2002
#      Create LDE_Example2.S
#
# ---------------------------------------------------------------------

options(width=100)

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

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

# ---------------------------------------------------------------------
# Here we add some simulated measurement error.  
# You'll want to remove this code if you're analyzing your own data.
    
simError1 <- rnorm(length(simData),mean=0,sd=0.5)
simError2 <- rnorm(length(simData),mean=0,sd=0.5)
simError3 <- rnorm(length(simData),mean=0,sd=0.5)
simError4 <- rnorm(length(simData),mean=0,sd=0.5)
L1 <- 1
L2 <- .5
L3 <- .4
L4 <- .7
uniqueX1 <- 01.25
uniqueX2 <- 01.75
uniqueX3 <- 01.0
uniqueX4 <- 02.25

theData <- cbind((L1 * simData) + simError1,
                 (L2 * simData) + simError2,
                 (L3 * simData) + simError3,
                 (L4 * simData) + simError4)

# ---------------------------------------------------------------------
# Here you set up the parameters for the analysis.

theInterval <- 1   # Time scale (interval between occasions)
theTimes <- 5      # Number of occasions used for estimation
theVars <- 4       # Number of indicator variables
minTau <- 1        # Minimum analysis lag
maxTau <- 16       # Maximum analysis lag

# ---------------------------------------------------------------------
# Here variables are initialized and a matrix is created to store the results.

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

resultsMatrixLDE1 <- matrix(NA, length(theTaus), 13)

system("rm resultsLDE1.dat")
system("rm temp.out")


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

tIndex <- 1

for (theTau in theTaus) {
    tFile <- paste("temp2.mx", sep="")
    cat("#define nocc ", theTimes, "   ! number of occasions of measurement\n", 
        file=tFile, append=F, sep="")
    cat("#define nind ", theVars, "   ! number of indicators\n", 
        file=tFile, append=T, sep="")
    cat("#define ni ", theTimes*theVars, "    ! Product nocc*nind\n", 
        file=tFile, append=T, sep="")
    cat("#define no ", totalN, "   ! sample size\n", 
        file=tFile, append=T, sep="")
    cat("#define time ", theInterval*theTau, "   ! measurement interval * number of lags\n", 
        file=tFile, append=T, sep="")

    tMatrix <- matrix(NA, theTimes, 3)
    for (i in 1:theTimes) {
        tMatrix[i,1] <- 1
        tMatrix[i,2] <- i - ((theTimes + 1)/2)
        tMatrix[i,3] <-(tMatrix[i,2])^2 / 2
    }
    tFile <- paste("temp4.mx", sep="")
    cat("  Matrix T time\n", file=tFile, append=F, sep="")
    cat("  Matrix K\n", file=tFile, append=T, sep="")
    write(t(format(round(tMatrix,2))), file=tFile, ncol=3, append=T)

    system("cat rungeAuto1.mx temp2.mx rungeAuto3.mx temp4.mx rungeAuto5.mx > temp.mx")

    tx <- matrix(NA, totalN, theTimes*theVars)
    for (i in 1:(totalN-((maxTimes-1)*theTau))) {
        for (k in 0:(maxTimes-1)) {
            tx[i,(0*maxTimes)+k+1] <- theData[(i+(k*theTau)),1]
            tx[i,(1*maxTimes)+k+1] <- theData[(i+(k*theTau)),2]
            tx[i,(2*maxTimes)+k+1] <- theData[(i+(k*theTau)),3]
            tx[i,(3*maxTimes)+k+1] <- theData[(i+(k*theTau)),4]
        }
    }
    tCov <- cov(tx, use="complete.obs")
    tFile <- paste("temp.cov", sep="")
    cat("*\n", file=tFile, append=F, sep="")
    write(t(format(round(tCov,5))), file=tFile, ncol=theTimes*theVars, append=T)
    system("mx < temp.mx >> temp.out") 
    tPars <- scan("t1pars.dat", 
            list(eta=0,zeta=0,covxdx=0,varx=0,vardx=0,SDd2x=0,chi2=0,
                b1=0,b2=0,b3=0,b4=0,b5=0,b6=0), skip=1)

    cat("\n Tau=", theTau, "\n")
    cat(cbind(tPars$chi2, tPars$zeta, tPars$eta, 
        tPars$varx, tPars$vardx, tPars$SDd2x, tPars$covxdx, 
        tPars$b1, tPars$b2, tPars$b3, tPars$b4, tPars$b5), "\n")
    
    t1 <- cbind(cbind(tPars$chi2, tPars$zeta, tPars$eta, 
        tPars$varx, tPars$vardx, tPars$SDd2x, tPars$covxdx, 
        tPars$b1, tPars$b2, tPars$b3, tPars$b4, tPars$b5))
    resultsMatrixLDE1[tIndex, 2:13] <- t1
    resultsMatrixLDE1[tIndex, 1] <- theTau
    write(t(format(round(resultsMatrixLDE1[tIndex,],5))), file="resultsLDE1.dat",
            ncol=13, append=T)
    tIndex <- tIndex + 1
}


colNames <- c("Tau", "ChiSq", "EstZeta", "EstEta", "VarX", "VarDX", "VarD2X", "CovXdX", "L1", "L2", "L3", "L4", "L5")
rowNames <- rep(NULL, dim(resultsMatrixLDE1)[1])
dimnames(resultsMatrixLDE1) <- list(rowNames,colNames)

round(resultsMatrixLDE1, 3)

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

# q()