R/MARSS_marxss.R

In gragusa/MARSS: Multivariate Autoregressive State-Space Modeling

###################################################################################
# marxss form definition file
# MARSS.marxss takes MARSS() call input 
# and returns the MARSS.call list with 2 model objects added
#    marssMODEL(form=marss) in the  $marss element
#    marssMODEL(form=marxss) in the $model element
# marss_to_marxss converts marssMODEL(form=marss) to marssMODEL(form=marxss)
# marxss_to_marss converts marssMODEL(form=marxss) to marssMODEL(form=marss)
# plus the following methods helper functions
# print_marxss, coef_marxss, predict_marxss, describe_marxss, MARSSinits_marxss
###################################################################################

###################################################################################
# Coerce model list input in MARSS() call to marxss model object (class marxss) put in $alt.forms
# Coerce marxss model into marss model and put in $marss
# using form = MARXSS
# x(t)=B(t)x(t-1) + U(t) + C(t)c(t) + w(t), W~MVN(0,Q)
# y(t)=Z(t)x(t) + A(t) + D(t)d(t) + v(t), V~MVN(0,R)
# x(t0) = x0 + l, L ~ MVN(0,V0)
# produces model object with fixed, free, tinitx, diffuse and data
# rhs is specified by fixed and free; lhs is specified with data (really should be list(x=, y=))
# attributes of model object has model.dims, X.names, form, equation

# The conversion functions has 3 parts
# Part 1 Set up the defaults and allowed structures
# Part 2 Do the conversion of model list to a marxss object
# Part 3 Do the conversion of marxss object to marss object

###################################################################################
MARSS.marxss=function(MARSS.call){
  #Part 1 Set up defaults and check that what the user passed in is allowed
  
  #Check that no args were passed into MARSS that are not allowed
  marxss.allowed.in.MARSS.call=c("model")
  allowed.in.call=c(marxss.allowed.in.MARSS.call,common.allowed.in.MARSS.call)
  if(any(!(names(MARSS.call) %in% allowed.in.call))){
    bad.names=names(MARSS.call)[!(names(MARSS.call) %in% allowed.in.call)]
    msg=paste("Argument ", paste(bad.names, collapse=", "),"  not allowed MARSS call for form ", MARSS.call$form, ". See ?MARSS.marxss\n", sep="")
    cat("\n","Errors were caught in MARSS.marxss \n", msg, sep="")
    stop("Stopped in MARSS.marxss() due to problem(s) with model specification.\n", call.=FALSE)
  }
  
  # 1 Check for form dependent user inputs for method and reset defaults for inits, MCbounds, and control if desired
  
  # 2 Specify the text shortcuts and whether factors or matrices can be passed in
  #   The names in the allowed list do not need to be A, B, Q .... as used in form=marss object
  #   Other names can be used if you want the user to use those names; then in the MARSS.form function
  #   you convert the user passed in names into the form marss names with the A, B, Q, R, ... names
  #   checkModelList() will check what the user passes in against these allowed values, so
  #   so you need to make sure each name in model.defaults has a model.allowed value here 
  model.allowed = list(
    A=c("scaling", "unconstrained", "unequal", "equal", "zero"),
    D=c("unconstrained", "equal", "unequal", "zero","diagonal and unequal","diagonal and equal","zero"),
    B=c("identity", "zero", "unconstrained", "unequal", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
    Q=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
    R=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
    U=c("unconstrained", "equal", "unequal", "zero"),
    C=c("unconstrained", "equal", "unequal", "zero","diagonal and unequal","diagonal and equal","zero"),
    x0=c("unconstrained", "equal", "unequal", "zero","diagonal and unequal","diagonal and equal"),
    V0=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
    Z=c("identity","unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "onestate"),
    c=c("zero"),
    d=c("zero"),
    tinitx=c(0,1),
    diffuse=c(TRUE,FALSE),
    factors = c("Z"),
    matrices = c("A","B","Q","R","U","x0","Z","V0","D","C","d","c")
  )
  
  #model.defaults is form dependent so you must specify it
  model.defaults =list(Z="identity", A="scaling", R="diagonal and equal", B="identity", U="unconstrained", Q="diagonal and unequal", x0="unconstrained", V0="zero", D="zero", d=matrix(0,1,1), C="zero", c=matrix(0,1,1), tinitx=0, diffuse=FALSE)
  
  if(!is.null(MARSS.call[["model"]][["c"]]))
    if(!identical(MARSS.call$model$c, "zero") & !all(MARSS.call$model$c==0)) model.defaults$C="unconstrained"
  if(!is.null(MARSS.call[["model"]][["d"]]))
    if(!identical(MARSS.call$model$d, "zero") & !all(MARSS.call$model$d==0)) model.defaults$D="unconstrained"
  
  #This checks that what user passed in model list can be interpreted and converted to form marss
  #if no errors, it updates the model list by filling in missing elements with the defaults
  MARSS.call$model=checkModelList( MARSS.call$model, model.defaults, model.allowed )
  
  # Part 2 Convert the model list to a marssMODEL object, form=marxss
  ## set up fixed and free elements
  fixed = free = list()
  
  model = MARSS.call[["model"]]
  model.elem = c("Z","A","R","B","U","Q","x0","V0","D","d","C","c") 
  dat = MARSS.call[["data"]]
  if(is.vector(dat)) dat=matrix(dat,1,length(dat))
  n = dim(dat)[1]; TT = dim(dat)[2]
  if(is.null(rownames(dat))){
    Y.names = paste("Y",seq(1, n),sep="") #paste(seq(1, n), sep="") 
    rownames(dat)=Y.names
  }else{ 
    Y.names = rownames(dat)
    if(any(duplicated(Y.names))){
      for(i in Y.names[duplicated(Y.names)]){
        loc = (Y.names==i)
        nn = sum(loc)
        Y.names[loc]=paste(i,"-",1:nn,sep="")
        rownames(dat)[loc] = Y.names[loc]
        MARSS.call[["data"]] = dat
      }
    }
  }
  
  ## Set m based on Z specification IF Z was specified; errors will be reported later if m conflicts with other parameters
  m = NA
  if (identical(model$Z, "unconstrained")) m = n
  if (identical(model$Z, "equalvarcov")) m = n
  if (identical(model$Z, "diagonal and equal")) m = n
  if (identical(model$Z, "diagonal and unequal")) m = n
  if (identical(model$Z, "onestate")) m = 1
  if (identical(model$Z, "identity")) m = n
  if (is.factor(model$Z)) m = length(levels(model$Z)) 
  if (is.array(model$Z)) m = dim(model$Z)[2] 
  
  X.names=NULL
  if(!(is.null(model[["X.names"]]))) X.names = model[["X.names"]]
    if(is.null(X.names) & identical(model$Z,"identity"))
      X.names=paste("X.",Y.names,sep="")
  if( is.null(X.names) & is.array(model$Z)){
    if(length(dim(model$Z))==3)
      if(dim(model$Z)[3]==1)
        if(is.design(model$Z) & !is.null(colnames(model$Z)))
          X.names=colnames(model$Z)
  }
  if( is.null(X.names) & is.factor(model$Z)){
    X.names=unique(as.character(model$Z))
  }
  if(is.null(X.names) & is.matrix(model$Z))
    if(is.design(model$Z) & !is.null(colnames(model$Z)))X.names=colnames(model$Z)
  if(is.null(X.names) & is.matrix(model$Z))
    if(is.identity(model$Z)) 
      X.names = paste("X.",Y.names,sep="")
  if(is.null(X.names)) X.names = paste("X",seq(1, m),sep="") #paste(seq(1, m),sep="")  #
  
  ## Set c1 based on C specification if C specified with a particular shape
  ## error checking later will complain if C and c (or D and d) conflict
  if (is.array(model$C)) c1 = dim(model$C)[2] else c1 = 1
  if (is.array(model$D)) d1 = dim(model$D)[2] else d1 = 1
  
  for(el in c("c","d")){
    thedim=get(paste(el,"1",sep=""))
    if(identical(model[[el]], "zero")){ model[[el]]=matrix(0,thedim,1) }
    #if(!any(is.na(model[[el]])) & all(model[[el]]==0)) model[[toupper(el)]]="zero"
    if(is.vector(model[[el]])) model[[el]]=matrix(model[[el]],1,length(model[[el]]))
  }
  #Now set c1 and d1 based on c and d, which should now be a matrix of some sort.  This ensures that c1 and d1 are set
  #if c and C or d and D conflict, this will be picked up the error checking because dims of C or D and model.dims won't match
  c1=dim(model$c)[1]; d1=dim(model$d)[1]
  #3rd dim of params are set to 1 and will be reset to correct value at end
  model.dims = list(data=c(n,TT),x=c(m,TT),y=c(n,TT),w=c(m,TT),v=c(n,TT),Z=c(n,m,1),U=c(m,1,1),A=c(n,1,1),B=c(m,m,1),Q=c(m,m,1),R=c(n,n,1),x0=c(m,1,1),V0=c(m,m,1),D=c(n,d1,1), d=c(d1,1,1), C=c(m,c1,1), c=c(c1,1,1))
  
  ## Error-checking section that is specific to marxss form
  # Note most error checking happens in checkMARSSInputs, checkModelList, and is.marssMLE
  #If a model elements passed in as a factors, make sure it is the correct length otherwise construction of marssMODEL object will break
  problem = FALSE
  msg=NULL
  ## Check model structures that are passed in as factors
  correct.factor.len = list(Z=n,A=n,R=n,B=m,U=m,Q=m,x0=m,V0=m)  
  for (el in model.elem) {
    #if a factor then it needs to have the correct length otherwise construction of marssMODEL object will break and no NAs
    if( is.factor(model[[el]]) ) {
      if(length(model[[el]]) != correct.factor.len[[el]]) {
        problem=TRUE
        msg = c(msg, paste(" The model$", el, " is being passed as a factor and should be length ", correct.factor.len[[el]], " based on data dims. It's not. See help file.\n", sep="")) }
      if(NA.in.fac <- NA %in% model[[el]]) {
        problem=TRUE
        msg = c(msg, paste(" NAs are not allowed in model factor for ", el, ". See help file.\n", sep="")) }
    }  #is factor  
  } # end for (el in model.elem)
  
  #if el == Z then factor needs to have m levels
  if( is.factor(model$Z) ) {
    if(length(levels(model$Z)) != m) {
      problem=TRUE
      msg=c(msg," When Z is a factor, the number of levels must equal the number of state processes (m).\n")
    } }
  
  #Check that if A is scaling, then Z spec must lead to a design matrix
  if(identical(model$A,"scaling")){
    if(is.array(model$Z) & length(dim(model$Z))==3)
      if(dim(model$Z)[3]!=1 && !is.design(model$Z, zero.cols.ok = TRUE)) { #if it is a matrix
        problem=TRUE
        msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\nYou can construct a scaling A matrix and pass that in.\n")
      }
    if((!is.array(model$Z) & !is.factor(model$Z))) #if it is a string
      if(!(model$Z %in% c("onestate","identity"))){
        problem=TRUE
        msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\n")
      }
    if(is.matrix(model$Z) && !is.design(model$Z, zero.cols.ok = TRUE)) { #if it is a matrix, won't be array due to first test
      problem=TRUE
      msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\n")
    }
  }  
  if(is.array(model$x0) & length(dim(model$x0))==3)
    if(dim(model$x0)[3]!=1){
      problem=TRUE
      msg = c(msg, " x0 cannot be time-varying thus if x0 in model arg is 3D, the 3rd dim must equal 1.\n")
    } 
  if(is.array(model$V0) & length(dim(model$V0))==3)
    if(dim(model$V0)[3]!=1){
      problem=TRUE
      msg = c(msg, " V0 cannot be time-varying thus if V0 in model arg is 3D, the 3rd dim must equal 1.\n")
    } 
  #if C is diagonal and equal or diagonal and unequal, then d1=m
  if(identical(model$C,"diagonal and equal") | identical(model$C,"diagonal and unequal"))
    if(c1!=m){
      problem=TRUE
      msg = c(msg, " If C is diagonal, it must be square and c must be m x 1.\n")
    }
  #if D is diagonal and equal or diagonal and unequal, then d1=n
  if(identical(model$D,"diagonal and equal") | identical(model$D,"diagonal and unequal"))
    if(d1!=n){
      problem=TRUE
      msg = c(msg, " If D is diagonal, it must be square and d must be n x 1.\n")
    }
  #if c and d can't have any NAs or Infs
  for(el in c("c","d"))
    if(any(is.na(model[[el]])) | !is.numeric(model[[el]]) | any(is.infinite(model[[el]]))){
      problem=TRUE
      msg = c(msg, paste(" ",el,"must be numeric and have no NAs, NaNs, or Infs.\n"))
    }
  #c and d must be a 2D matrix and 2nd dim must be 1 or TT
  for(el in c("c","d")){
    if( length(dim(model[[el]]))!=2 ){
      problem=TRUE
      msg = c(msg, paste(" ",el,"must be a 2D matrix.\n"))
    }else{
      if( !(dim(model[[el]])[2] == 1 | dim(model[[el]])[2] == TT) ){
        problem=TRUE
        msg = c(msg, paste(" ",el,"must be a 2D matrix with 2nd dim equal to 1 or T (length of data).\n"))
      }
    }
  }
  
  #If there are problems
  if(problem)  {
    cat("\n","Errors were caught in MARSS.marxss \n", msg, sep="")
    stop("Stopped in MARSS.marxss() due to problem(s) with model specification.\n", call.=FALSE)
  }
  #end of error section  
  
  #Change c and d to array so that it can be handled by the normal fixed/free constructions
  for(el in c("c","d")){
    row.names=rownames(model[[el]])
    model[[el]]=array(model[[el]],dim=c(dim(model[[el]])[1],1,dim(model[[el]])[2]))
    rownames(model[[el]])=row.names
  }
  
  
  ##Translate the text shortcuts into a marssMODEL object 
  ## Translate the model structure names (shortcuts) into fixed and free
  ## fixed is a dim(1)*dim(2) X 1 vector of the fixed (intercepts) values
  ## free is a dim(1)*dim(2) X p vector of the free (betas) values for the p estimated elements
  
  model.elem = c("Z","A","R","B","U","Q","x0","V0","D","C","d","c")
  if(which(model.elem=="Z")>which(model.elem=="A")) model.elem=rev(model.elem) #Z must go first
  
  tmp=list()
  for(el in model.elem) {
    if(MARSS.call$silent==2) cat(paste("  Building fixed and free matrices for ",el,".\n",sep=""))
    tmp[[el]]="not assigned"
    if(el=="Z" & is.factor(model$Z)) { 
      tmp[[el]] = matrix(0,model.dims$Z[1], model.dims$Z[2])  
      for(i in X.names) tmp[[el]][which(model$Z==i), which(as.vector(X.names)==i)] = 1
    }
    if(el=="Z" & identical(model$Z,"onestate")) {   #m=1
      tmp[[el]] = matrix(1, n, 0)
    }
    if( identical(model[[el]],"identity") ) { 
      tmp[[el]] = diag(1,model.dims[[el]][1])
    }
    if(identical(model[[el]],"diagonal and equal")) {
      tmp[[el]] = array(list(0),dim=c(model.dims[[el]][1],model.dims[[el]][2]))
      diag(tmp[[el]])="diag" #paste(el,"(diag)",sep="")
      if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
    }
    if(identical(model[[el]],"diagonal and unequal")) {
      tmp[[el]] = array(list(0),dim=c(model.dims[[el]][1],model.dims[[el]][2]))
      dim.mat = model.dims[[el]][1]
      el.labs=as.character(1:dim.mat)
      if(el %in% c("V0","Q","B")) el.labs=X.names
      if(el %in% c("Z","R")) el.labs=Y.names
      diag(tmp[[el]])=paste("(",el.labs,",", el.labs,")",sep="") #paste(el,"(",as.character(1:dim.mat),",",as.character(1:dim.mat),")",sep="")
      if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
    }
    if(identical(model[[el]],"unconstrained") | identical(model[[el]],"unequal")){
      tmp[[el]]=array(NA,dim=c(model.dims[[el]][1],model.dims[[el]][2])) 
      if(el %in% c("Q","R","V0")){  #variance-covariance matrices
        dim.mat = model.dims[[el]][1]
#         for(i in 1:dim.mat){
#           for(j in 1:dim.mat) tmp[[el]][i,j]=tmp[[el]][j,i]=paste("(",i,",",j,")",sep="") #paste(el,"(",i,",",j,")",sep="")
#         }
          tmp[[el]]=matrix(paste("(",rep(1:dim.mat,dim.mat),",",rep(1:dim.mat,each=dim.mat),")",sep=""),dim.mat,dim.mat)
          tmp[[el]][upper.tri(tmp[[el]])] = t(tmp[[el]])[upper.tri(t(tmp[[el]]))]  
      }else{ #not var-cov matrix
        row.name=1:model.dims[[el]][1]
        col.name=1:model.dims[[el]][2]
        if(el %in% c("U","x0")) row.name=X.names
        if(el == "A") row.name=Y.names
        if(el %in% c("C","D")){
          if(el=="C") row.name=X.names
          if(el=="D") row.name=Y.names
          if(!is.null(rownames(model[[tolower(el)]]))) col.name=rownames(model[[tolower(el)]])
        }
#         for(i in 1:model.dims[[el]][1]){
#           for(j in 1:model.dims[[el]][2]){
#             if(model.dims[[el]][2]>1) tmp[[el]][i,j]=paste("(",row.name[i],",",col.name[j],")",sep="") #paste(el,"(",row.name[i],",",col.name[j],")",sep="")
#             else tmp[[el]][i,j]=paste(row.name[i],sep=",") #paste(el,row.name[i],sep=",")
#           }
#         }
        if(model.dims[[el]][2]>1){
          tmp[[el]]=matrix(paste("(",rep(row.name,model.dims[[el]][2]),",",rep(col.name,each=model.dims[[el]][1]),")",sep=""),model.dims[[el]][1],model.dims[[el]][2])
        }else{
          tmp[[el]]=matrix(row.name,model.dims[[el]][1],1)
        }
      }
      if(length(tmp[[el]])==1) tmp[[el]][1,1]=el 
    } #unconstrained
    if(identical(model[[el]],"equalvarcov")) {
      tmp[[el]]=array("offdiag",dim=c(model.dims[[el]][1],model.dims[[el]][2])) #array(paste(el,"(offdiag)",sep=""),dim=model.dims[[el]])
      diag(tmp[[el]])="diag" #paste(el,"(diag)",sep="")
      if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
    }
    if(identical(model[[el]],"equal")) { 
      tmp[[el]]=array("1",dim=c(model.dims[[el]][1],model.dims[[el]][2])) #array(el,dim=model.dims[[el]])
    }
    if(identical(model[[el]],"zero")) { 
      tmp[[el]]=array(0,dim=c(model.dims[[el]][1],model.dims[[el]][2]))
    }
    if(is.array(model[[el]])) {
      tmp[[el]]=model[[el]]
    }
    if(el=="A" & identical(model[[el]], "scaling")) {  #check above ensures that Z is design and time-invariant
      ## Construct A from fixed Z matrix
      tmp[[el]] = matrix(list(),model.dims$A[1],model.dims$A[2])
      tmp[[el]][,1]=Y.names  
      for(i in 1:m) {
        nonzeroZ=tmp$Z[,i]!=0
        if(any(nonzeroZ)) tmp[[el]][min(which(nonzeroZ)), 1] = 0
      }
    }     
    if(identical(tmp[[el]],"not assigned")) stop(paste("MARSS.marxss: tmp was not assigned for ",el,".\n",sep=""))
    tmpconst=convert.model.mat(tmp[[el]])
    free[[el]] = tmpconst$free
    fixed[[el]] = tmpconst$fixed 
    
    #set the last dim of the model.dims since it was at a temp value to start
    model.dims[[el]][3]=max(dim(free[[el]])[3],dim(fixed[[el]])[3])
  }
  
  #save the row names for the inputs by setting in fixed
  for(el in c("c","d")){
    if(is.null(rownames(tmp[[el]]))) rownames(tmp[[el]])=paste(el,seq(1,dim(tmp[[el]])[1]),sep="")
    rownames(fixed[[el]])=rownames(tmp[[el]])
  }
  
  
  #Set the marssMODEL form marxss
  #This is the f+Dp form for the MARXSS model used for user displays, printing and such
  marxss_object = list(fixed=fixed, free=free, data=dat, tinitx=model$tinitx, diffuse=model$diffuse)
  #set the attributes
  class(marxss_object) = "marssMODEL"
  attr(marxss_object, "obj.elements")=c("fixed","free","data","tinitx","diffuse")
  attr(marxss_object, "form")="marxss"
  attr(marxss_object, "model.dims")=model.dims
  #par.names are what needs to be in fixed/free pair
  attr(marxss_object, "par.names")=c("Z","A","R","B","U","Q","x0","V0","D","C","d","c")
  attr(marxss_object, "X.names")=X.names
  attr(marxss_object, "Y.names")=Y.names
  attr(marxss_object, "equation")="x_{t}=B_{t}*x_{t-1}+U_{t}+C_{t}*c_{t}+w_{t}; w_{t}~MVN(0,Q_{t})\ny_{t}=Z_{t}*x_{t}+A_{t}+D_{t}*d_{t}+v_{t}; v_{t}~MVN(0,R_{t})"

  #Change alldefaults global to match the form
  alldefaults[[MARSS.call$method]][["inits"]][["C"]]=0
  alldefaults[[MARSS.call$method]][["inits"]][["D"]]=0
  assign("alldefaults", alldefaults, pkg_globals)
  
  ## Check that the marssMODEL object output by MARSS.form() is ok since marxss_to_marss will go south otherwise
  if(MARSS.call$silent==2) cat(paste("  Running is.marssMODEL on the marxss model.\n",sep=""))
  tmp = is.marssMODEL(marxss_object, method=MARSS.call$method)
  if(!isTRUE(tmp)) {
    cat(tmp) 
    stop("Stopped in MARSS.marxss() due to problem(s) with model specification.", call.=FALSE)
  }
  
  #Put the marxss model into $model since model holds the model in the 'form' form
  MARSS.call$model = marxss_object
  
  ## Create marssMODEL(form=marss) object added to call
  #when called with a marssMODEL object (as here), marxss_to_marss returns a marssMODEL object
  MARSS.call$marss=marxss_to_marss(marxss_object)
  
  ## Return MARSS call list with $marss and $model added
  MARSS.call
}

marss_to_marxss=function(x, C.and.D.are.zero=FALSE){
  if(!(class(x) %in% c("marssMODEL","marssMLE"))) stop("marss_to_marxss: this function needs a marssMODEL or marssMLE object")
  #This function returns a MLE object where the model and par parts of the MLE object are in marxss form for printing purposes.
  #This function needs a marxss marssMODEL object and will break otherwise
  #You cannot back construct a marxss from the marss model
  #The function will also work is x is a model object, but then it just returns marxss.marssMODEL
  #written this way so it doesn't crash if x happens to be a marssMODEL in case I later dynamically write function names/calls
  
  if(class(x)=="marssMODEL"){
    marss.model=x
    if(!("marss" %in% attr(marss.model,"form"))) stop("marss_to_marxss: this function requires a marssMODEL object in marss form.\n",call.=FALSE)
  }else{ 
    marss.model=x[["marss"]]
    if(!("marss" %in% attr(marss.model,"form"))) stop("marss_to_marxss: this function requires a marssMLE object with element $marss a marssMODEL in marss form.\n",call.=FALSE)
  }
  
  if(!C.and.D.are.zero){
    if(class(x)=="marssMODEL"){ 
      stop("marss_to_marxss: function was called with a marss model object instead of MLE object, so needs a marxss model passed in.")
    }else{ 
      marxss.model=x[["model"]] #should be model since we want the marxss form
      if(any(is.null(attr(marxss.model,"form")),!("marxss" %in% attr(marxss.model,"form")))){
        stop("marss_to_marxss: function was called with a MLE object, so needs the $model element to be form marxss.")
      }
    }
  }else{ #C and D are zero so we can construct a marxss object
    marxss.model = marss.model #use the marss model as a template
    marxss.dims=attr(marss.model,"model.dims")
    marxss.model[["fixed"]][["C"]]=array(0,dim=c(marxss.dims[["x"]][1],1,1))
    marxss.model[["fixed"]][["D"]]=array(0,dim=c(marxss.dims[["y"]][1],1,1))
    marxss.model[["free"]][["C"]]=array(0,dim=c(marxss.dims[["x"]][1],0,1))
    marxss.model[["free"]][["D"]]=array(0,dim=c(marxss.dims[["y"]][1],0,1))
    marxss.model[["fixed"]][["c"]]=array(0,dim=c(1,1,1))
    marxss.model[["fixed"]][["d"]]=array(0,dim=c(1,1,1))
    marxss.model[["free"]][["c"]]=array(0,dim=c(1,0,1))
    marxss.model[["free"]][["d"]]=array(0,dim=c(1,0,1))
    marxss.dims[["C"]]=c(marxss.dims[["x"]][1],1,1)
    marxss.dims[["D"]]=c(marxss.dims[["y"]][1],1,1)
    marxss.dims[["c"]]=c(1,1,1)
    marxss.dims[["d"]]=c(1,1,1)
    #reset the attributes to marxss form
    #obj.elements, X.names and Y.names stay the same
    attr(marxss.model, "form")="marxss"
    attr(marxss.model, "model.dims")=marxss.dims
    #par.names are what needs to be in fixed/free pair; order is important
    attr(marxss.model, "par.names")=c("Z","A","R","B","U","Q","x0","V0","D","C","d","c")
    attr(marxss.model, "equation")="x_{t}=B_{t}*x_{t-1}+U_{t}+C_{t}*c_{t}+w_{t}; w_{t}~MVN(0,Q_{t})\ny_{t}=Z_{t}*x_{t}+A_{t}+D_{t}*d_{t}+v_{t}; v_{t}~MVN(0,R_{t})"
  }
  if(class(x)=="marssMODEL") return(marxss.model) #in marxss form

  x[["model"]]=marxss.model #now in marxss form
  marxss.dims = attr(marxss.model, "model.dims")
  
  #got here, so class of x is marssMLE
    for(val in c("par","start","par.se","par.bias","par.upCI","par.lowCI")){
      if(!is.null(x[[val]])){
          tmp.dim=dim(marxss.model$free$C)[2] #how many C parameters
        if(tmp.dim==0){
          x[[val]][["C"]] = matrix(0,0,1)
        }else{
          #because marss.U is [marxss.C marxss.U]
          x[[val]][["C"]] = x[[val]][["U"]][1:tmp.dim,, drop=FALSE]
          x[[val]][["U"]] = x[[val]][["U"]][-(1:tmp.dim),, drop=FALSE]
        }
          tmp.dim=dim(marxss.model$free$D)[2] #how many D parameters
        if(tmp.dim==0){
          x[[val]][["D"]] = matrix(0,0,1)
        }else{
          #because marss.A is [marxss.D marxss.A]
          x[[val]][["D"]] = x[[val]][["A"]][1:tmp.dim,, drop=FALSE]
          x[[val]][["A"]] = x[[val]][["A"]][-(1:tmp.dim),, drop=FALSE]
        }
        for(el in c("c","d")){
          x[[val]][[el]] = matrix(0,0,1) #because c and d are inputs not estimated
        }
      } #not is null val
    } #for val in par, start

    #returning a marssMLE object where the par, start etc are in marxss form
  # $model is in marxss and $marss stays the same
  return(x) 
}

marxss_to_marss=function(x, only.par=FALSE){
  #x is a marssMODEL of form marxss
  #this will create a marss model object (if !only.par) and a par list in form marss
  #if only.par=TRUE then only the par element is changed and marss is used for the marss object
  
  #hold on to this since x will be changing and need to know what to return
  class.x=class(x)
  if(!(class.x %in% c("marssMODEL","marssMLE"))){
    stop("marxss_to_marss: x$model must be a marssMODEL or marssMLE.",call.=FALSE)
  }
  
  #check form if user passed in marssMODEL
  if(class.x=="marssMODEL"){
    if(!("marxss" %in% attr(x, "form"))) stop("marxss_to_marss: this function requires a marssMODEL object in marxss form.")
  }
  
  if(class.x=="marssMLE"){ #Then set the par elements to correspond to marss if they are in marxss form
    #check that the model element they passed in is marxss
    if(!("marxss" %in% attr(x[["model"]], "form"))) stop("marxss_to_marss: x$model must be in marxss form.",call.=FALSE)
    
    x.marss = list()
    
    #this will go through x and reset any par-like obj in marxss form
    for(val in c("par","start","par.se","par.bias","par.upCI","par.lowCI")){
      if(!is.null(x[[val]])){
          #because marss.U is [marxss.C marxss.U]
          x.marss[[val]][["U"]] = rbind(x[[val]][["C"]],x[[val]][["U"]])
          #because marss.A is [marxss.D marxss.A]
          x.marss[[val]][["A"]] = rbind(x[[val]][["D"]],x[[val]][["A"]])
          #other elements are the same as for marxss
          for(el in c("R","Q","B","Z","x0","V0"))
            x.marss[[val]][[el]]=x[[val]][[el]]
          #reset x[[val]] so it only includes the marss elements
          x[[val]]=x.marss[[val]] #replace the x[[val]] with marss version
      } #not is null val
    } #for val in par, start
    
    marxss.model=x[["model"]]
  }else{ #end if x is marssMLE
    marxss.model = x
  }
  #marxss.model is a marssMODEL in marxss form and x is the original marssMLE object
  #if only.par updating was requested, return x
  if(class.x == "marssMLE" & only.par & !is.null(x[["marss"]])) return(x)
  
  #else construct a marss model object from marxss.model and put in $marss
  
  marxss.dims=attr(marxss.model,"model.dims")  #fixed and free will be modified, so these are holders not shortcuts
  fixed=marxss.model[["fixed"]]
  free=marxss.model[["free"]]
  #marss.dims will be modified, so this is a holder not a shortcut
  marss.dims=marxss.dims
  n=marss.dims[["y"]][1]; m=marss.dims[["x"]][1]; TT=marss.dims[["x"]][2]
  #This step converts U+Cc into equivalent Uu and A+Dd into Aa
  #So U --> [C U] and u --> [c \\ 1]; A --> [D A] and a --> [d \\ 1]
  #This code adds fixed$u and fixed$a, and changes fixed$U and fixed$A; otherwise all stays the same
  for(el in c("C","D")){
    if(el=="C") el2="U" else el2="A"
    #if el all zero (fixed and all zero), it doesn't appear in the equation and U-->U and u-->1
    if(!is.fixed(marxss.model[["free"]][[el]]) | !all(sapply(marxss.model[["fixed"]][[el]],function(x){isTRUE(all.equal(x,0))}))){      
      #create fixed$u or fixed$a by add 1 to bottom of fixed$c or fixed$d
      #since fixed$c is a 3D array, we need to do this in an odd way:
      fixed[[tolower(el2)]]=array(apply(fixed[[tolower(el)]],3,rbind,1),dim=dim(fixed[[tolower(el)]])+c(1,0,0))
      #this fixed$u is just an input so we set the free to not estimated
      free[[tolower(el2)]]=array(0,dim=c(1,0,1)) #not estimated so 0 columns
      
      #next change U to [C U] and A to [D A]; need to figure out the dims since
      #C or U might be time-varying
      Tmax.fixed=max(dim(fixed[[el]])[3],dim(fixed[[el2]])[3])
      Tmax.free=max(dim(free[[el]])[3],dim(free[[el2]])[3])
      #dim. is c(dim 1 of C, dim 2 of C, dim 1 of U, dim 2 of U)
      dim.fixed=c(dim(fixed[[el]])[1],dim(fixed[[el]])[2],dim(fixed[[el2]])[1],dim(fixed[[el2]])[2])
      dim.free=c(dim(free[[el]])[1],dim(free[[el]])[2],dim(free[[el2]])[1],dim(free[[el2]])[2])
      #now that the dimensions are known, create and array holder for fixed$U and fixed$A
      tmp.fixed=array(NA,dim=c(dim.fixed[1]+dim.fixed[3],1,Tmax.fixed))
      tmp.free=array(0,dim=c(dim.free[1]+dim.free[3],dim.free[2]+dim.free[4],Tmax.free))
      #the first rows of fixed$U are fixed$C
      tmp.fixed[1:dim.fixed[1],,]=fixed[[el]]
      #the next rows are marxss.model$fixed$U
      tmp.fixed[(dim.fixed[1]+1):dim(tmp.fixed)[1],,]=fixed[[el2]]
      #Now create the new free$U.  If C estimated, free$C appears in the upper left
      if(dim.free[2]>0) tmp.free[1:dim.free[1],1:dim.free[2],]=free[[el]]
      #If U estimated, free$U appears in the lower right
      if(dim.free[4]>0) tmp.free[(dim.free[1]+1):dim(tmp.free)[1],(dim.free[2]+1):dim(tmp.free)[2],]=free[[el2]]
      #retain the column names (estimated parameter names)
      colnames(tmp.free)=c(colnames(free[[el]]),colnames(free[[el2]]))
      #assign the new fixed$U and free$U to fixed and free
      fixed[[el2]]=tmp.fixed
      free[[el2]]=tmp.free
      #settign U and A dims
      marss.dims[[el2]][2]=marxss.dims[[el]][2]+marxss.dims[[el2]][2]
    }else{  #Both C and U are all zero (fixed and all zero) 
      #so u is just 1 (a 1x1 matrix)
      fixed[[tolower(el2)]]=array(1,dim=c(1,1,1))
      free[[tolower(el2)]]=array(0,dim=c(1,0,1)) #not estimated
    }
  }
  
  #Now the fixed/free specify x=Bx+U(t)u(t)+w(t) and y=Zx+A(t)a(t)+v(t)
  #this part converts U(t)u(t) to U(t) and A(t)a(t) to A(t)
  #This requires making a vec(U(t)) that is specified by (rhs at end):
  #vec(Uu)=(t.u kron I_m)vec(U)=fixed+free*p=(t.u kron I_m)f+(t.u kron I_m)Dp
  #where f and D are fixed$U and free$U for U(t)u(t) form
  #the small case are inputs and the large case are estimated parameters  
  for(el in c("U","A")){
    #if c or a passed in. If not will be array(1,dim=c(1,1,1))
    #this if statement is just avoiding unneccesary code.  The math should still hold whether or not c is 1
    if(!identical(unname(fixed[[tolower(el)]]), array(1,dim=c(1,1,1)))){
      #hold onto fixed$U and free$U (not marxss.model$fixed and free but the new ones)
      free.orig=free[[el]]; fixed.orig=fixed[[el]]
      dim.free2=dim(free.orig)[2]; dim.free3=dim(free.orig)[3];  dim.fixed3=dim(fixed.orig)[3]; 
      dim.u.3=dim(fixed[[tolower(el)]])[3]
      Tmax=max(dim.free3, dim.fixed3, dim.u.3)
      #need the new marss U and A dims here which were defined above
      free[[el]]=array(0,dim=c(marss.dims[[el]][1],dim.free2,Tmax))
      colnames(free[[el]])=colnames(free.orig)
      fixed[[el]]=array(0,dim=c(marss.dims[[el]][1],1,Tmax))
      for(t in 1:Tmax){
        #the f and D of U (or A)
        f.t=sub3D(fixed.orig,t=min(t,dim.fixed3))
        d.t=sub3D(free.orig,t=min(t,dim.free3))
        #column vector of the u (or a) at time t
        ua.t=sub3D(fixed[[tolower(el)]],t=min(t,dim.u.3))  
        #vec(Uu)=(t.u kron I_m)vec(U)=(t.u kron I_m)f+(t.u kron I_m)Dp
        #again we want the new marss.dims defined above
        free[[el]][,,t]=(t(ua.t)%x%diag(1,marss.dims[[el]][1]))%*%d.t
        fixed[[el]][,,t]=(t(ua.t)%x%diag(1,marss.dims[[el]][1]))%*%f.t
      }
    }
  }
  marss.elem = c("Z","A","R","B","U","Q","x0","V0")
  free=free[marss.elem]
  fixed=fixed[marss.elem]
  dim3s=apply(rbind(unlist(lapply(free[marss.elem],function(x){dim(x)[3]})), unlist(lapply(fixed[marss.elem],function(x){dim(x)[3]}))),2,max)
  marss.dims = list(data=c(n,TT),x=c(m,TT),y=c(n,TT),w=c(m,TT),v=c(n,TT),
                    Z=c(n,m,dim3s[["Z"]]),U=c(m,1,dim3s[["U"]]),A=c(n,1,dim3s[["A"]]),
                    B=c(m,m,dim3s[["B"]]),Q=c(m,m,dim3s[["Q"]]),R=c(n,n,dim3s[["R"]]),
                    x0=c(m,1,1),V0=c(m,m,1))
  
  ## Create the marss marssMODEL object
  marss.model = list(fixed=fixed, free=free, data=marxss.model[["data"]], tinitx=marxss.model[["tinitx"]], diffuse=marxss.model[["diffuse"]])
  #set the attributes that change
  class(marss.model) = "marssMODEL"
  attr(marss.model, "form")="marss"
  attr(marss.model, "obj.elements")=c("fixed","free","data","tinitx","diffuse")
  attr(marss.model, "model.dims")=marss.dims
  attr(marss.model, "par.names")=marss.elem
  attr(marss.model, "X.names")=attr(marxss.model,"X.names")
  attr(marss.model, "Y.names")=attr(marxss.model,"Y.names")
  attr(marss.model, "equation")="x_{t}=B_{t}*x_{t-1}+U_{t}+w_{t}; w_{t}~MVN(0,Q_{t})\ny_{t}=Z_{t}*x_{t}+A_{t}+v_{t}; v_{t}~MVN(0,R_{t})"
  if(class.x=="marssMODEL"){ return(marss.model) #marssMODEL of form marss
  }else{ 
    #class.x=marssMLE, then adds the marss element to the marssMLE object
    x[["marss"]]=marss.model
    return(x) #returning a marssMLE object where the par, start etc are in marss form and marss element is set
  }
}

#the par element of a marssMLE object is in form=marss.  Convert to form=marxss for printing
print_marxss = function(x){ return(marss_to_marxss(x)) } 

#the par element of a marssMLE object is in form=marss.  Convert to form=marxss for printing
coef_marxss = function(x){ return(marss_to_marxss(x)) } #this uses $model for marxss object

MARSSinits_marxss = function(MLEobj, inits){
  if(is.null(MLEobj[["model"]])){
    stop("MARSSinits_marxss: this function needs a marssMODEL in marxss form in $model",call.=FALSE)
  }else{
    if(class(MLEobj[["model"]])!="marssMODEL") stop("MARSSinits_marxss: this function needs a marssMODEL in marxss form in $model",call.=FALSE)
    if(!("marxss" %in% attr(MLEobj[["model"]],"form"))) stop("MARSSinits_marxss: this function needs a marssMODEL in marxss form in $model",call.=FALSE)    
  }
    
  #B, Z, R, Q, x0 and V0 stay the same
  #U and A change
  #this function will return a U and A element for inits
  if(is.null(inits)) inits=list()
  inits.defaults=list(
    U=alldefaults[[MLEobj$method]][["inits"]][["U"]], C=0,
    A=alldefaults[[MLEobj$method]][["inits"]][["A"]], D=0 )
  elems=c("U","A","C","D")
  for(elem in elems){
    tmp.dim=dim(MLEobj$model$free[[elem]])[2] #how many estimated pars in marxss vers
    if(!is.null(inits[[elem]])){
      if(!(length(inits[[elem]]) %in% c(tmp.dim,1))) 
        stop(paste("MARSSinits: ", elem," inits must be either a scalar (dim=NULL) or a matrix with 1 col and rows equal to the num of est values in ",elem,".",sep=""), call.=FALSE )
      if(tmp.dim!=0) inits[[elem]] = matrix(inits[[elem]],tmp.dim,1) else inits[[elem]]=matrix(0,0,1)
    }else{  inits[[elem]] = matrix(inits.defaults[[elem]],tmp.dim,1) }
  }
  inits$U = rbind(inits$C,inits$U) #yes, C on top
  inits$A = rbind(inits$D,inits$A) #yes, D on top
  return(inits)
}

predict_marxss = function(x, newdata, n.ahead, t.start){
  #x is a marssMLE object
  
  #This takes the newdata argument from a predict call and interprets the inputs in the context of the form of x$model
  #It will return a marssMODEL (form=marss) object ready for use in prediction
  #n.ahead is 1 by default; t.start is TT+1 by default
  #this makes tinitx=0, and uses E(x)_(t.start-1) as if t.start=1, init x is specified by x0T; V0 is V0T
  #if t.start!=1, init x is specified by xtT[,t.start-1] and VtT[,,t.start-1]
  
  #First get a par list in marxss form
  marxss.par = marss_to_marxss(x)[["par"]] #we only need par changed since marxss is in $model
  marxss.dims = attr(x[["model"]], "model.dims") #the marxss model dims
  marxss.model = x[["model"]]
  if( !("marxss" %in% attr(marxss.model,"form")) ) stop("predict_marxss: x$model needs to be in marxss form.", call.=FALSE)
  TT=marxss.dims[["y"]][2]
  form="marxss"
  allow.in=c("data","c","d") #allowed in newdata
  
  if(!(class(newdata) %in% c("list","data.frame")))
    stop("predict_marxss: newdata must be a list or dataframe.",call.=FALSE)
  
  #next interpret newdata;
  
  #if user passes in a data.frame, make an effort to interpret that
  #the following code makes sure newdata is a list with data, c and d with same 3rd dim (n.ahead)
  if(is.data.frame(newdata)){ #try to construct the list
    newdata.dataframe=newdata
    newdata=list()
    names.dataframe=names(newdata.dataframe)
    if(any(duplicated(names.dataframe))) stop("predict_marxss: the dataframe should not have any duplicated names",call.=FALSE)
    
    #first construct the data matrix from the dataframe
    Y.names=attr(marxss.model,"Y.names")
    #find any column names in the dataframe that match the rownames in the data matrix
    Y.match=match(Y.names,names.dataframe)
    if(any(!is.na(Y.match))){
      if(dim(newdata)[1]!=n.ahead){
        stop("predict_marxss: you have passed data in with newdata as a dataframe.\nThe number of rows must equal n.ahead in this case.",call.=FALSE)
      }
      cat(paste("Alert: y (data) are present in the dataframe and prediction will be conditioned on these values.\n",collapse=""))
      newdata[["data"]]=newdata.dataframe[Y.names[!is.na(Y.match)]]
    }
    #those that are missing are replaced with NA
    if(any(is.na(Y.match)))
      newdata[["data"]][Y.names[is.na(Y.match)]]=as.numeric(NA)
    #make sure the columns are in the same order as Y.names
    newdata[["data"]]=newdata$data[Y.names]
    newdata[["data"]]=t(as.matrix(newdata[["data"]])) #time across columns
    
    if(!(dim(newdata)[1]==n.ahead | dim(newdata)[1]!=1)){
      stop("predict_marxss: The number of rows in the newdata dataframe must be 1 or n.ahead.",call.=FALSE)
    }
    
    #Create the c and d matrices
    for(el in c("c","d")){
      is.zero = all(marxss.model[["fixed"]][[toupper(el)]]==0) & all(marxss.par[[toupper(el)]]==0)
      if(is.zero){ #if C or D is all zero then c or d is not needed; make it all 0
        newdata[[el]]=matrix(0,marxss.dims[[el]][1],1)
        rownames(newdata[[el]])=rownames(marxss.model[["fixed"]][[el]])
      }else{ #need c or d
        el.names=rownames(marxss.model[["fixed"]][[el]])        
        #find any column names in the dataframe that match the rownames in the el
        el.match=match(el.names,names.dataframe)
        #subset those columns in the dataframe that match the el names
        el.in.dataframe=el.names[!is.na(el.match)]
        newdata[[el]]=newdata.dataframe[el.in.dataframe]
        #make it into a matrix with time going across the columns like in $model
        newdata[[el]]=t(as.matrix(newdata[[el]]))
        #deal with any missing c or d rows
        if(!all(el.names %in% names.dataframe)){
          bad.names = el.names[!(el.names %in% names.dataframe)]
          if((n.ahead+t.start-1)>TT){
            #require that user passes in all the c and d inputs                       
            stop(c("predict_marxss: some of the ",el," inputs are missing: ",paste(bad.names,collapse=", ")),call.=FALSE)
          }else{ #replace missing names with values from the model
            if(marxss.dims[[el]][3]==1) t.el=1 else t.el=t.start:(t.start+n.ahead-1)
            tmp.el=matrix(marxss.model[["fixed"]][[el]][bad.names,1,t.el],length(bad.names),t.start+n.ahead-1)
            rownames(tmp.el)=bad.names
            newdata[[el]]=rbind(newdata[[el]],tmp.el)
          }
        }      
        #makesure the columns are in the same order as el.names
        newdata[[el]]=newdata[[el]][el.names,,drop=FALSE]      
        if(any(is.na(newdata[[el]]))){
          stop(paste("predict_marxss: there cannot be any NAs in the ",el," matrix.\n",sep=""),call.=FALSE)
        }
      }
    } #for el
  } #newdata is now a list with elements data, c, and d
  
  #if user passed in a list
  if(is.list(newdata)){
    for(el in c("c","d")){
      is.zero = all(marxss.model[["fixed"]][[toupper(el)]]==0) & all(marxss.par[[toupper(el)]]==0)
      if(is.zero){ #if C or D is all zero then c or d is not needed; make it all 0
        newdata[[el]]=matrix(0,marxss.dims[[el]][1],1)
        rownames(newdata[[el]])=rownames(marxss.model[["fixed"]][[el]])
      }else{ #need c or d
        if(!(el %in% names(newdata))){
          if((n.ahead+t.start-1)>TT){
            #require that user passes in all the c and d inputs                       
            stop(c("predict_marxss: Model has ", toupper(el), " but ", el," is missing from newdata\n and cannot be inferred since prediction extends beyond original dataset."),call.=FALSE)
          }else{ #replace missing names with values from the model
            if(marxss.dims[[el]][3]==1) t.el=1 else t.el=t.start:(t.start+n.ahead-1)
            newdata[[el]]=matrix(marxss.model[["fixed"]][[el]][,,t.el],marxss.dims[[el]][1],t.start+n.ahead-1)
            rownames(newdata[[el]])=rownames(marxss.model[["fixed"]][[el]])
          }
        }
        if(!is.matrix(newdata[[el]]))
          stop(c("predict_marxss: ", el," in newdata must be a matrix with time across the columns and n.ahead columns."),call.=FALSE)
        names.list=rownames(newdata[[el]])
        el.names=rownames(marxss.model[["fixed"]][[el]])
        #find any row names in the matrix that match the rownames in the el
        el.match=match(el.names,names.list)
        #subset those columns in the dataframe that match the el names
        el.in.matrix=el.names[!is.na(el.match)]
        newdata[[el]]=newdata[[el]][el.in.matrix,,drop=FALSE]
        if(!all(el.names %in% names.list)){
          bad.names = el.names[!(el.names %in% names.list)]
          if((n.ahead+t.start-1)>TT){
            #require that user passes in all the c and d inputs                       
            stop(c("predict_marxss: some of the ",el," inputs are missing: ",paste(bad.names,collapse=", ")," \nand cannot be inferred since prediction is beyond the end of the original dataset."),call.=FALSE)
          }else{ #replace missing names with values from the model
            if(marxss.dims[[el]][3]==1) t.el=1 else t.el=t.start:(t.start+n.ahead-1)
            tmp.el=matrix(marxss.model[["fixed"]][[el]][bad.names,1,t.el],length(bad.names),t.start+n.ahead-1)
            rownames(tmp.el)=bad.names
            newdata[[el]]=rbind(newdata[[el]],tmp.el)
          }
        }      
        #makesure the columns are in the same order as el.names
        newdata[[el]]=newdata[[el]][el.names,,drop=FALSE]      
        if(any(is.na(newdata[[el]]))){
          stop(paste("predict_marxss: There cannot be any NAs in the ",el," matrix.\n",sep=""),call.=FALSE)
        }
      }
    }
    if(("data" %in% names(newdata))){
      el="data"
      Y.names=attr(marxss.model,"Y.names")
      el.names=Y.names
      if(!is.matrix(newdata[[el]]))
        stop(c("predict_marxss: ", el," in newdata must be a matrix with time across the columns and n.ahead columns."),call.=FALSE)
      
      names.list=rownames(newdata[[el]])
      if(!all(el.names %in% names.list)){
        bad.names = el.names[!(el.names %in% names.list)]
        stop(c("predict_marxss: Some of the ",el," rows are missing: ",paste(bad.names,collapse=", ")),call.=FALSE)
      }
      #find any row names in the matrix that match the rownames in the el
      el.match=match(el.names,names.list)
      #subset those columns in the dataframe that match the el names
      el.in.matrix=el.names[!is.na(el.match)]
      newdata[[el]]=newdata[[el]][el.in.matrix,]
      #makesure the columns are in the same order as el.names
      newdata[[el]]=newdata[[el]][el.names,] 
    }
  }
  
  
  #now do some error checking
  if(!all(names(newdata) %in% allow.in)) stop(paste("predict_marxss: only allowed inputs for form ",form," are ",allow.in,".",sep=""),call.=FALSE)
  #Check that the dims of inputs are the same (or = 1) and fix c or d that have length 1 to have same dim 2 as other inputs
  thedims=unlist(lapply(newdata,function(x){dim(x)[2]}))
  #only consider those with length!=1
  thedims=thedims[thedims!=1]
  #n.ahead will have some value, since n.ahead is 1 by default in predict_marssMLE which called this function.
  #make sure n.ahead and any dim 2 of data, c or d !=1 are the same
  thedims=c(n.ahead,thedims)
  if(!all( abs(thedims - mean(thedims)) < .Machine$double.eps )){
    stop(paste("predict_marxss: in newdata, the 2nd dimension of all inputs (and n.ahead if passed in) must be equal (if not =1).",sep=""))
  }
  
  #set data to be all missing if it wasn't passed in
  if(is.null(newdata[["data"]])){
    newdata[["data"]] = matrix(as.numeric(NA), marxss.dims[["y"]][1], n.ahead)
    rownames(newdata[["data"]]) = attr(marxss.model,"Y.names")
  }else{ #it was passed in
    if(dim(newdata[["data"]])[2]!=n.ahead)
      stop(c("predict_marxss: data in newdata must be a matrix with time across the columns and n.ahead columns."),call.=FALSE)
  }
  
  #Now newdata ready as a list with elements c and d that correspond to model c and d as passed into MARSS()
  
  #### make a list of time-varying parameters
  param.names = attr(marxss.model,"par.names")
  time.varying.fixed = c()
  time.varying.free = c()
  for( elem in param.names ) {
    if( dim(marxss.model[["free"]][[elem]])[3] == 1 ){
      time.varying.free[elem] = FALSE  #not time-varying
    }else{ time.varying.free[elem] = TRUE }
    if( dim(marxss.model[["fixed"]][[elem]])[3] == 1 ){
      time.varying.fixed[elem] = FALSE  #not time-varying
    }else{ time.varying.fixed[elem] = TRUE }
  }
  time.varying = time.varying.fixed | time.varying.free
  
  if(any(time.varying)){
    if((t.start+n.ahead-1)>TT)
      stop(paste("predict_marxss: ",paste(param.names[time.varying], collapse=", ")," are time-varying.\nIn this case, you cannot forecast past the end of the time series\n(t.start+n.ahead must be < length of original data).\n",sep=""),call.=FALSE)
    param.t=t.start:(t.start+n.ahead-1)
  }
  #Now we need to construct a marssMODEL object (form=marxss) for predicting
  pred.marxss=marxss.model #copy the elements and attributes that won't change
  pred.model.dims=marxss.dims
  pred.model.dims[["data"]][2]=n.ahead
  pred.model.dims[["x"]][2]=n.ahead
  pred.model.dims[["y"]][2]=n.ahead
  for( el in param.names[!(param.names %in% c("c","d"))] ){
    if(!time.varying.fixed[el]){
      pred.marxss[["fixed"]][[el]]=marxss.model[["fixed"]][[el]]
    }else{ #is time.varying
      pred.marxss[["fixed"]][[el]]=marxss.model[["fixed"]][[el]][,,param.t,drop=FALSE]
      pred.model.dims[[el]][3]=n.ahead
    }
    if(!time.varying.free[el]){
      pred.marxss[["free"]][[el]]=marxss.model[["free"]][[el]]
    }else{ #is time.varying
      pred.marxss[["free"]][[el]]=marxss.model[["free"]][[el]][,,param.t,drop=FALSE]
      pred.model.dims[[el]][3]=n.ahead
    }  }
  
  for(el in c("c","d")){
    #this forces 2nd dim of matrix to be n.ahead
    pred.model.dims[[el]][3]=n.ahead
    #now replace any c or d row passed in with the new values
    pred.marxss[["fixed"]][[el]] = array(newdata[[el]],dim=c(dim(newdata[[el]])[1],1,n.ahead))
    rownames(pred.marxss[["fixed"]][[el]]) = rownames(newdata[[el]])
    #by definition free c and d is not estimated
    pred.marxss[["free"]][[el]] = marxss.model[["free"]][[el]]
  }
  
  #Set the initial conditons
  #Set these to the estimated distribution of x at t.start-1 conditioned on all the data (original)
  pred.marxss[["tinitx"]]=0
  kf=MARSSkf(x) #x is the original marssMLE obj passed in
  if(t.start==1){
    pred.marxss$fixed$x0=array(kf$x0T, dim=marxss.dims$x0)
    pred.marxss$fixed$V0=array(vec(kf$V0T), dim=c(marxss.dims$V0[1]*marxss.dims$V0[2],1,1))
  }else{
    pred.marxss$fixed$x0=array(kf$xtT[,t.start-1], dim=marxss.dims$x0)
    pred.marxss$fixed$V0=array(vec(kf$VtT[,,t=(t.start-1)]), dim=c(marxss.dims$V0[1]*marxss.dims$V0[2],1,1))
  }
  pred.marxss$free$x0=array(0,dim=c(marxss.dims$x0[1],0,1))
  pred.marxss$free$V0=array(0,dim=c(marxss.dims$V0[1]*marxss.dims$V0[2],0,1))
  
  pred.marxss[["data"]]=newdata[["data"]]
  attr(pred.marxss, "model.dims")=pred.model.dims #the marxss model dims
  
  ##Next step is to create the marssMLE object ready for predict.  This is done by changing x

  #These are not estimated when x_(tstart-1) and V_(t.start-1) are fixed at expected values
  for(val in c("par","start")){
    for(el in c("x0","V0")){
      if(!is.null(x[[val]][[el]])){
        x[[val]][[el]]=matrix(0,0,1)
      }
    }
  }
  #These don't have meaning when x_(tstart-1) and V_(t.start-1) are fixed at expected values
  for(val in c("par.se","par.bias","par.upCI","par.lowCI")){
    for(el in c("x0","V0")){
      if(!is.null(x[[val]][[el]])){
        x[[val]][[el]]=NULL
      }
    }
  }
  
  
  #now pred.marxss is a marssMODEL object (form=marxss) to be used for prediction
  x[["model"]]=pred.marxss
  #convert pred.marxss to a marss object and put in $marss
  x[["marss"]]=marxss_to_marss(pred.marxss)
  #only pass in marssMODEL since x par has been corrected above
  
  #x is now a marssMLE object with its marss and model changed to reflect newdata
  # with y, c and d potentially changed and t.start and t.end changed and TT changed
  # par$x0 and $V0 changed to reflect t.start
 
  #return this object to predict.marssMLE; only MLEobj is needed
  #but newdata list returned for debugging
  return(list(MLEobj=x, newdata=newdata))
}

describe_marxss = function(modelObj){ describe_marss(modelObj) }
#describe_marss works generally with marxss form models (of which marss is one type)
#describe_marss is in the file describe_marssMODEL.R


########################################################################
# is.marssMODEL_marxss function
# Check that the marxss object has all the parts it needs
# fixed, free, and par.names
# and that these have the proper size and form
# m is pulled from fixed$x0
########################################################################
is.marssMODEL_marxss <- function(modelObj, method="kem"){
  
  msg=NULL

  if( !("marxss" %in% attr(modelObj, "form") ) ){ 
    msg = c(msg, "Form attribute of the model object does not include marxss.\n")
  }
  ## Set up par.names that should be marxss model
  en = c("Z","A","R","B","U","Q","x0","V0","D","C","d","c")
  
  #Check that par.names has these and only these names
  par.names=attr(modelObj, "par.names")
  if( !all(en %in% par.names ) ){ 
    msg = c(msg, "Element ", en[!(en %in% par.names)], " is missing from the par.names attribute of the model object.\n")
  }
  if( !all( par.names %in% en ) ) { 
    msg = c(msg, "Only ", en, "should be in the par.names attribute of the model object.\n")
  }
  if(!is.null(msg)){  #rest of the tests won't work so stop now
    return(msg)
  }
  
  ###########################
  # Check model.dims are correct
  ###########################
  n = dim(modelObj$data)[1]
  TT = dim(modelObj$data)[2]
  m = dim(modelObj$fixed$x0)[1]
  c1=dim(modelObj$fixed$c)[1]
  d1=dim(modelObj$fixed$d)[1]
  en = c("Z","A","R","B","U","Q","x0","V0","D","C","d","c", "data", "x", "y", "w", "v")
  correct.dim1 = c(Z=n,A=n,R=n,B=m, U=m, Q=m, x0=m, V0=m, D=n, C=m, c=c1, d=d1, data=n, x=m, y=n, w=m, v=n)
  correct.dim2 = c(Z=m,A=1,R=n,B=m, U=1, Q=m, x0=1, V0=m, D=d1,C=c1, c=TT, d=TT, data=TT, x=TT, y=TT, w=TT, v=TT)
  model.dims=attr(modelObj, "model.dims")
  for (elem in en) {
    ## Check for problems in the fixed/free pairs. Problems show up as TRUE 
    dim.flag1 = dim.flag2 = FALSE
    
    # check dim
    dim.flag1 = c(dim.flag1,!( model.dims[[elem]][1] == correct.dim1[[elem]] )) 
    dim.flag2 = c(dim.flag2,!( model.dims[[elem]][2] == correct.dim2[[elem]] ))
  }  
  if (any(c(dim.flag1, dim.flag2))) {  #There's a problem
    if(any(dim.flag1)) {
      msg = c(msg, paste("Dim 1 of ", en[dim.flag1], "is incorrect. Dims should be ", correct.dim1[dim.flag1], ", for a marss model.\n"))
    }
    if(any(dim.flag2)) {
      msg = c(msg, paste("Dim 2 of ", en[dim.flag2], "is incorrect. Dims should be ", correct.dim2[dim.flag2], ", for a marss model.\n"))
    }
    msg=c("\nErrors were caught in is.marssMODEL_marxss()\n", msg)
    return(msg)
  }  
  
  fixed=modelObj$fixed; free=modelObj$free
  
  ###########################
  # Check that x0 and V0 are not time-varying
  ###########################
  en = c("x0", "V0")
  time.var = NULL
  for (elem in en) {
    time.var.flag = FALSE
    time.var.flag = dim(fixed[[elem]])[3]!=1 | dim(free[[elem]])[3]!=1
    time.var <- c(time.var, time.var.flag)
  }
  if(any(time.var)) {  #There's a problem
    msg = c(msg, paste(en[time.var], "cannot be time-varying.  3rd dim of fixed and free must equal 1.\n"))
    msg=c("\nErrors were caught in is.marssMODEL_marxss()\n", msg)
    return(msg)
  }
  
  ###########################
  # Check that none of the var-cov matrices have negative values on the diagonal
  # and that there are no f+Dq elements only f+0q or 0+Dq
  # and D must be a design matrix, so no beta_1*q1 + beta_2*q2 elements
  ###########################
  en = c("R", "Q", "V0")
  neg = bad.var = not.design = NULL
  for (elem in en) {
    neg.flag = bad.var.flag = not.design.flag = FALSE
    for(i in 1:max(dim(free[[elem]])[3],dim(fixed[[elem]])[3])){
      if(dim(fixed[[elem]])[3]==1){i1=1}else{i1=i}
      if(dim(free[[elem]])[3]==1){i2=1}else{i2=i}
      if(is.fixed(free[[elem]][,,min(i,dim(free[[elem]])[3]),drop=FALSE])){ #this works on 3D mats
        zero.free.rows = matrix(TRUE,correct.dim1[[elem]]*correct.dim2[[elem]],1)
      }else{
        zero.free.rows=apply(free[[elem]][,,i2,drop=FALSE]==0,1,all) #works on 3D mat
        #the requirement is that each estimated element (in p) appears only in one place in the varcov mat, but fixed rows (0 rows) are ok
        not.design.flag = !is.design(free[[elem]][,,i2,drop=FALSE],strict=FALSE,zero.rows.ok=TRUE,zero.cols.ok=TRUE) #works on 3D if dim3=1
      }
      zero.fixed.rows=apply(fixed[[elem]][,,i1,drop=FALSE]==0,1,all) #works on 3D
      fixed.mat = unvec(fixed[[elem]][,,i1],dim=c(correct.dim1[[elem]],correct.dim2[[elem]]))
      if( any(!zero.fixed.rows & !zero.free.rows) ){
        bad.var.flag = TRUE   #no f+Dq rows
      }
      if(any(takediag(fixed.mat)<0,na.rm=TRUE)) neg.flag=TRUE      #no negative diagonals
    } #end the for loop over time
    not.design = c(not.design, not.design.flag)
    neg = c(neg, neg.flag)
    bad.var = c(bad.var, bad.var.flag)
  } #enf the for loop over elem
  if(any(neg)) {
    msg = c(msg, paste("Negative values are on the diagonal of ", en[neg], ". Neg values are illegal on the diag of a var-cov matrix.\n", sep=""))
  }
  if(any(bad.var)) {
    msg = c(msg, paste("Fixed and estimated values are combined in some elements of ", en[bad.var], ". This is not allowed.\n", sep=""))
  }
  if(any(not.design)) {
    msg = c(msg, paste("The D matrices of ", en[not.design], " must be design matrices.\n", sep=""))
  }
  
  ###########################
  # Check that V0, Q and R matrices are symmetric and positive-definite
  ###########################
  en = c("R", "Q", "V0")
  pos = symm = NULL
  for (elem in en) {
    varcov.flag = TRUE; varcov.msg=""
    var.dim = c(correct.dim1[[elem]],correct.dim2[[elem]])
    for(i in 1:model.dims[[elem]][3]){
      if(dim(fixed[[elem]])[3]==1){i1=1}else{i1=i}
      if(dim(free[[elem]])[3]==1){i2=1}else{i2=i}
      #works on 3D if dim3=1
      par.as.list = fixed.free.to.formula(fixed[[elem]][,,i1,drop=FALSE],free[[elem]][,,i2,drop=FALSE],var.dim) #coverts the fixed,free pair to a list matrix
      tmp=is.validvarcov(par.as.list, method=method)
      varcov.flag=varcov.flag & tmp$ok
      if(!tmp$ok) varcov.msg = c(varcov.msg, paste(" ", tmp$error, "at t=", i, "\n",sep=""))
      
      if(!varcov.flag) msg = c(msg, paste("The variance-covariance matrix ", elem, " is not properly constrained.\n", sep=""), varcov.msg)
    } #end for loop over time

  } #end for loop over elements

if(length(msg) == 0){ return(NULL)
}else {
  msg=c("\nErrors were caught in is.marssMODEL_marxss()\n", msg)
  return(msg)
}
}