# ---------------------------------------------------------------------
# Program: LDE_2ndOrder_OneIndicator_5D.R
#  Author: Steven Boker and Pascal Deboeck
#    Date: Sat Sep  5 12:23:32 EDT 2009
# This program runs an example single indicator 2nd order LDE model
#   with 5 time-delay columns in the embedded state space matrix
# ---------------------------------------------------------------------
# Revision History
#    -- Sat Sep  5 12:23:32 EDT 2009
#      Created untitled.
# 2018-12-09 Updated for modern fit function and new data URL
# ---------------------------------------------------------------------

# ---------------------------------------------------------------------
# Variables 
# ---------------------------------------------------------------------

# ----------------------------------
# Read libraries and set options.

require(OpenMx)

# ----------------------------------
# Read demo data.

# data(demoOscillator)

demoOscillator <- read.table("http://openmx.ssri.psu.edu/sites/default/files/demoOscillator.txt", col.names="x")

# ==================
# = Set constants. =
# ==================

tau    <- 1.0  # The lag between subsequent columns in the embedded matrix
deltaT <- 0.3  # The amount of time elapsed between subsequent observations
embedD <- 7.0  # The number of columns in the time-delay embedded matrix

# ===================================
# = Time delay embed the demo data. =
# ===================================

Embed <- function(x, E, tau) {
	# Create a time delay matrix from x with embedding dimension E and lag tau
	len <- length(x)
	out <- x[1:(len-(E*tau)+tau)]
	for(i in 2:E) {
		out <- cbind(out,x[(1+((i-1)*tau)):(len-(E*tau)+(i*tau))])
	}
	return(out)
}
	
embeddedOscillator <- Embed(demoOscillator$x, embedD, tau)
dimnames(embeddedOscillator) <- list(NULL, paste("x", 1:embedD, sep=""))

# ========================================
# = Create the fixed LDE loading matrix. =
# ========================================

L1 <- rep(1,embedD)
L2 <- c(1:embedD)*tau*deltaT-mean(c(1:embedD)*tau*deltaT)
L3 <- (L2^2)/2
LDE.Original <- cbind(L1,L2,L3)

# ===============================
# = Create 2nd order LDE model. =
# ===============================

manifestVars <- dimnames(embeddedOscillator)[[2]]

ldeModel1 <- mxModel("LDE_Model_1",
    mxMatrix(name="L", "Full", values=LDE.Original, free=FALSE, byrow=TRUE),
    mxMatrix(name="A", "Full", 3, 3, byrow=TRUE,
        values=c(  0,  0, 0,
                   0,  0, 0,
                 -.2,-.2, 0), 
        free=c(FALSE,FALSE,FALSE,
               FALSE,FALSE,FALSE,
                TRUE, TRUE,FALSE)
    ),
    mxMatrix(name="S", "Symm", 3, 3, byrow=TRUE,
        values=c(.8, -.1, 0,
                -.1,  .8, 0,
                  0,   0, .8), 
        free=c( TRUE, TRUE, FALSE,
               TRUE, TRUE,FALSE,
               FALSE,FALSE, TRUE), 
        lbound=c(0.000001, -10000  , 0.000001,
                 -10000  , 0.000001, 0.000001,
                 0.000001, 0.000001, 0.000001)
    ),
    mxMatrix(name="U", "Diag", embedD, embedD, values=.8,free=TRUE, lbound=0.000001),
    mxMatrix(name="I", "Iden", 3),
    mxAlgebra(name="R", L %*% solve(I-A) %*% S %*% t(solve(I-A)) %*% t(L) + U,
        dimnames = list(manifestVars, manifestVars)
    ),
    mxExpectationNormal("R"), mxFitFunctionML(),
    mxData(cov(embeddedOscillator), type="cov", numObs=dim(embeddedOscillator)[1])
)

# ======================================================
# = Fit the LDE model and examine the summary results. =
# ======================================================

ldeModel1Fit <- mxRun(ldeModel1)

summary(ldeModel1Fit)