R/dse1.R

In dse: Dynamic Systems Estimation (Time Series Package)

###
###  There is still some function description here and in EvalEst that
###    should be moved to help files.
###

.onAttach  <- function(libname, pkgname) {
   .DSEflags(list(COMPILED=TRUE))
   invisible(TRUE)
   }

.DSEflags <- local({
   .DSE.flags <- character(0)
   function(new) {
       if(!missing(new))
           .DSE.flags <<- new
       else
           .DSE.flags
       }
})

# With NAMESPACEs use .onAttach instead of .First.Lib 
#.First.lib <- function(library, section) {
#   assign(".DSEflags()$COMPILED", TRUE, env = .GlobalEnv)
#   assign(".DSEDUP",      TRUE, env = .GlobalEnv)
#   invisible(library.dynam("dse"))
#   }
  

DSEversion <- function() 
  {if (!is.R()) return("version cannot be determined.") else
   z <- c("setRNG", "tframe", "dse","EvalEst",  
          "curve", "CDNmoney", "tsfa", "TSdbi")
   z <- z[ z %in% library()$results[,1] ]
   r <- NULL
   for (pac in z ) 
     r <- c(r,  packageDescription(pac, fields="Version"))
   names(r) <- z
   r
  }


##############################################################################

# The code was divided roughly into 
#   the groups listed below, but the organization has changed a little bit.
#   The code grouping can be seen by grep ing on the string '<<<<<<' eg:
#    grep "<<<<<<" dse.R

#    Functions which work on a model (i.e. if a model with data is allowed as
#             an arguement then the data is ignored):
#        -model summary, description, print, comparison and 
#              calculation of properties
#        -model conversion 
             
#    Functions which work on data ( and a model ):
#        -likelihood and residual calculation
#        -statistical tests
#        -parameter estimation 
#        -model reduction (this could work just on models but
#                    comparisons require data)     

#    Utility Functions:
#        -utilities for generating theoretical statistics 
#              (i.e.- model generated, not from data)
#        -data transformations
#        -model and data scaling
#        -utilities for polynomial arithmetic
#        -internal utilities used for updating objects 
#        -data interface functions

#   Of special note in the internal utilities are two programs,
#   setArrays and setTSmodelParameters which take a model list and
#   update the representation and parameter information respectively. 
#   i.e. setArrays uses the parameter information and ensures that the 
#    array representation is consistent while setTSmodelParameters uses the 
#    arrays and ensures that the parameter vector is consistent.

#############################################################################
#############################################################################
residuals.TSestModel <- function(object, ...) object$estimates$pred-outputData(object)

# I don't think this exists really works for namespaces
 if (!exists("acf.default", mode="function")){
   acf.default  <- stats::acf
   acf <- function(x, ...)UseMethod("acf")
   }
 
 acf.TSdata <- function(x, ...)acf(outputData(x))
 acf.TSestModel <- function(x, ...) acf(x$estimates$pred - outputData(x), ...)

#############################################################################

#functions which work on models   <<<<<<<<<<


############################################################

#     functions for model summary, description, print and
#      comparison and functions for calculation of model properties

############################################################


print.TSestModel <- function(x, ...) 
{ cat("neg. log likelihood=",x$estimates$like[1],"\n")
  if(!is.null(x$converged)) {if(!x$converged) cat(" (not converged)") }
  print(x$model,...) 
  invisible(x)
}

print.SS <- function(x, digits=options()$digits, latex=FALSE, ...) 
    {#  (... further arguments, currently disregarded)
     printM <- function(st, M, digits, latex)
	{arrayc <- function(n)
	   {cat("\\left[ \\begin{array}{"); for(j in 1:n) cat("c"); cat("}\n") }
	 if (latex) cat("\\begin{equation}\n")
	 cat("\n", st, "=\n")
	 col <- if (is.matrix(M)) dim(M)[2] else length(M)
	 if (latex) 
	   {arrayc(col)
	    if (!is.matrix(M)) M <- matrix(M, 1, length(M))
	    M <- signif(M, digits=digits)
            for (i in 1:dim(M)[1])
	      {for (j in 1:dim(M)[2]) 
	         {cat(M[i,j])
		  if (j != dim(M)[2]) cat(" & ")
		 }
	       if (i != dim(M)[1]) cat("\\\\ \n")
	      }
	    cat("\n\\end{array} \\right] \n\\end{equation}\n")
	   }
	 else print(M, digits=digits)
	 invisible()
	}
	
     printM("F", x$F, digits, latex)
     if(!is.null(x$G)) printM("G", x$G, digits, latex)
     printM("H", x$H, digits, latex)
     if ( is.nonInnov.SS(x)) 
       {printM("Q", x$Q, digits, latex)
        printM("R", x$R, digits, latex)
       }
     else if (is.innov.SS(x))  printM("K", x$K, digits, latex)
     if(!is.null(x$z0)) printM("initial state", x$z0, digits, latex)
     if(!is.null(x$rootP0))
        printM("square root of initial state tracking error", x$rootP0, digits, latex)
     else if(!is.null(x$P0))
        printM("initial state tracking error", x$P0, digits, latex)
     invisible(x)
    }

print.ARMA <- function(x, digits=options()$digits, latex=FALSE, L=TRUE, fuzz=1e-10, ...) 
#  (... further arguments, currently disregarded)
   {arrayc <- function(n)
	{cat("\\left[ \\begin{array}{"); for(j in 1:n) cat("c"); cat("}\n") }
        
     A <- x$A
     B <- x$B
     C <- x$C
     if(!is.null(x$TREND))
        cat("TREND= ", format(signif(x$TREND,digits=digits)))
     if (latex | L)
       {if (latex) cat("\\begin{equation}")
        cat("\nA(L) =\n")
        if (latex) arrayc(dim(A)[3])
        for(i in 1:dim(A)[2]) 
          {for(j in 1:dim(A)[3]) 
             {cat(format(signif(A[1,i,j],digits=digits)))
              if(dim(A)[1] > 1) for(l in 2:dim(A)[1]) 
                 if (abs(A[l,i,j]) > fuzz)
                  {if(1==sign(A[l,i,j])) cat("+")
                   cat(format(signif(A[l,i,j],digits=digits)))
                   if (latex) cat("L^{",l-1,"}", sep="") else cat("L",l-1,sep="")
                  }
               if (latex & j != dim(A)[3]) cat(" & ") else cat("    ")
             }
           if (latex) cat("\\\\ \n") else cat("\n")
         }
        if (latex) cat("\\end{array} \\right] \n\\end{equation}\n\\begin{equation}")
        cat("\nB(L) =\n")
        if (latex) arrayc(dim(B)[3]) 
        for(i in 1:dim(B)[2]) 
          {for(j in 1:dim(B)[3]) 
             {cat(signif(B[1,i,j],digits=digits))
              if (2 <= dim(B)[1]) for(l in 2:dim(B)[1]) 
                 if (abs(B[l,i,j]) > fuzz)
                  {if(1==sign(B[l,i,j])) cat("+")
                   cat(signif(B[l,i,j],digits=digits))
                   if (latex) cat("L^{",l-1,"}", sep="") else cat("L",l-1,sep="")
                  }
              if (latex & j != dim(B)[3]) cat(" & ") else cat("    ")
             }
           if (latex) cat("\\\\ \n") else cat("\n")
         }
         if (latex) cat("\\end{array} \\right] \n\\end{equation}\n")
         if(!is.null(x$C)) 
          {if (latex) cat("\\begin{equation}")
           cat("\nC(L) =\n")
           if (latex) arrayc(dim(C)[3])
           for(i in 1:dim(C)[2]) 
           {for(j in 1:dim(C)[3]) 
             {cat(signif(C[1,i,j],digits=digits))
              if (2<=dim(C)[1]) for(l in 2:dim(C)[1]) 
                 if (abs(C[l,i,j]) > fuzz)
                  {if(1==sign(C[l,i,j])) cat("+")
                   cat(signif(C[l,i,j],digits=digits))
                   if (latex) cat("L^{",l-1,"}", sep="") else cat("L",l-1,sep="")
                  }
              if (latex & j != dim(C)[3]) cat(" & ") else cat("    ")
             }
             if (latex) cat("\\\\ \n") else cat("\n")
           } 
           if (latex) cat("\\end{array} \\right] \n\\end{equation}\n")
	  }
       }
     else
       {for(l in 1:dim(A)[1]) {cat("\nA(L=",l-1,")\n");print(A[l,,],digits=digits)}
        for(l in 1:dim(B)[1]) {cat("\nB(L=",l-1,")\n");print(B[l,,],digits=digits)}
        if(!is.null(x$C))
          for(l in 1:dim(C)[1]) {cat("\nC(L=",l-1,")\n");print(C[l,,],digits=digits)}
       }
     invisible(x)
} 
 
summary.TSestModel <- function(object, ...)
  {#  (... further arguments, currently disregarded)
   residual <- residuals(object)
   sampleT <- nrow(residual)
   p <- ncol(residual)	
   #Om <- t(residual) %*% residual/sampleT
   Om <- crossprod(residual) /sampleT
   rmse <- matrix( diag(Om)^.5 ,1,p)
   dimnames(rmse) <- list(c("RMSE"), seriesNamesOutput(object))

   classed(list(  # summary.TSestModel constructor
     estimates=list(
        l=object$estimates$like[1],
        rmse=rmse,
        sampleT=sampleT,
        converged= object$converged,
        nlmin.results= (!is.null(object$nlmin.results)),
        conv.type= object$nlmin.results$conv.type,
        filter= (!is.null(object$filter)),
        smooth= (!is.null(object$smooth)) ),
     model=summary(object$model)), "summary.TSestModel")
  }


print.summary.TSestModel <- function(x, digits=options()$digits, ...)
  {#  (... further arguments, currently disregarded)
   cat("neg. log likelihood =",x$estimates$l)
   cat("    sample length ="     ,x$estimates$sampleT, "\n")
   print(x$estimates$rmse)
   if (!is.null(x$estimates$converged))
                           cat("convergence: ",x$estimates$converged,"\n")
   if (x$estimates$nlmin.results)
                           cat("convergence type: ", x$estimates$conv.type,"\n")
   if (x$estimates$filter) cat("Includes  filter  estimates.\n")
   if (x$estimates$smooth) cat("Includes smoother estimates.\n")
   print(x$model, digits=digits)
   invisible(x)
  }



summary.SS <- function(object, ...)
  {#  (... further arguments, currently disregarded)
   m <- nseriesInput(object)
   p <- nseriesOutput(object)
   n <- nstates(object)
   classed(list(  # summary.SS constructor
         description=object$description,
         input.series=seriesNamesInput(object),
         output.series=seriesNamesOutput(object),
         innov=is.innov.SS(object),
         m=m,
         p=p,
         n=n,
         P=n * (m+2*p),  #assumes full rank noise
         P.actual = length(coef(object)),
         P.IC = sum(object$location %in% c("z", "P", "r P")),
	 constants=length(object$const),
         ICs=(!is.null(object$z0)),
         init.track=(!is.null(object$P0)) | (!is.null(object$rootP0))
	 ), "summary.SS")
  }

print.summary.SS <- function(x, digits=options()$digits, ...)
    {#  (... further arguments, currently disregarded)
     if (x$innov) cat("innovations form ")
     cat("state space: ")
     cat(x$description,"\n")
     cat("inputs : ", x$input.series, "\n")
     cat("outputs: ", x$output.series, "\n")
     cat("   input  dimension = ", x$m)
     cat("   state  dimension = ",x$n)
     cat("   output dimension = ",x$p,"\n")
     cat("   theoretical parameter space dimension = ",x$P,"\n")
     cat("  ",x$P.actual, " actual parameters")
     if (0 != x$P.IC) cat(" (of which ",x$P.IC, " are ICs)")
     cat("   ",x$constants," non-zero constants\n")
     if (x$ICs)        cat("   Initial values specified.\n")
     else              cat("   Initial values not specified.\n")
     if (!x$innov)
       {if (x$init.track) cat("   Initial tracking error specified.\n")
        else              cat("   Initial tracking error not specified.\n")
       }
     invisible(x)
    }

summary.ARMA <- function(object, ...)
  {#  (... further arguments, currently disregarded)
   m <- nseriesInput(object)
   p <- nseriesOutput(object)
   classed(list(  # summary.ARMA constructor
         description=object$description,
         input.series=seriesNamesInput(object),
         output.series=seriesNamesOutput(object),

         a=dim(object$A)[1]-1,
         b=dim(object$B)[1]-1,
         c=dim(object$C)[1]-1, 
         m=m,
         p=p,
         P.actual=length(coef(object)),
         constants=length(object$const),
         trend=(!is.null(object$TREND)) ), "summary.ARMA")
}
 
print.summary.ARMA <- function(x, digits=options()$digits, ...)
    {#  (... further arguments, currently disregarded)
     cat("ARMA: ")
     cat(x$description,"\n")
     cat("inputs : ", x$input.series, "\n")
     cat("outputs: ", x$output.series, "\n")
     cat("     input  dimension = ", x$m)
     cat("     output dimension = ", x$p,"\n")
     cat("     order A = ", x$a)
     cat("     order B = ", x$b)
     cat("     order C = ", x$c,"\n") 
     cat("     ",x$P.actual, " actual parameters")
     cat("     ",x$constants," non-zero constants\n")
     if(x$trend) cat("     trend estimated.\n")
     else        cat("     trend not estimated.\n")
     invisible(x)
    }
 
tfplot.TSestModel <- function(x, ..., 
    tf=NULL, start=tfstart(tf), end=tfend(tf), 
    select.inputs=NULL, select.outputs=NULL,
    Title=NULL, xlab=NULL, ylab=NULL, 
    graphs.per.page=5, mar=par()$mar, reset.screen=TRUE) {

  # plot one-step ahead estimates and actual data.
  # ... is a list of models of class TSestModel.
  model <- x
  if (is.null(Title))
     Title <- "One step ahead predictions (dotted) and actual data (solid)"
  p<- nseriesOutput(model)
  if (is.null(select.outputs)) select.outputs <-1:p
  if (all(0==select.outputs)) select.outputs <- NULL
  Ngraphs <- length(select.outputs)
  if (!is.null(select.inputs)) if (all(0==select.inputs)) select.inputs <- NULL
  m<- nseriesInput(model)
  if (is.null(m)) m <-0
  else Ngraphs <- Ngraphs+length(select.inputs)  # NULL is zero
  Ngraphs <- min(Ngraphs, graphs.per.page)
  if(reset.screen) {
     old.par <- par(mfcol = c(Ngraphs, 1), mar = mar, no.readonly=TRUE)
     on.exit(par(old.par)) }
  names <-seriesNamesOutput(model)
  if (m!=0) names <-c(seriesNamesInput(model), names)
  if (is.null(names)) names <- rep(" ", m+p)
  if (1 == length(xlab)) xlab <- rep(xlab, m+p)

  if (m!=0) 
    {for (i in select.inputs) 
      {z <-NULL 
       for (model in append(list(x),list(...)) )
           z<-tbind(z,inputData(model, series=i))
       tframe(z) <-tframe(inputData(model))
       tfOnePlot(z,start=start,end=end, xlab=xlab[i], ylab=names[i]) # tsplot
       if(!is.null(Title) && (i==1) && (is.null(options()$PlotTitles)
                || options()$PlotTitles)) title(main = Title)
    } }
  for (i in select.outputs ) 
    {#z <-c(outputData(model, series=i),
     #      rep(NA,Tobs(model$estimates$pred) - TobsOutput(model)))
     z <- outputData(model, series=i)
     for (model in append(list(x),list(...))){
         if (! is.TSestModel(model)) 
	    stop("Argument in ... to be plotted is not a TSestModel object.")
         #z <- cbind(z,model$estimates$pred[,i,drop=FALSE])}
         z <- tbind(z,selectSeries(model$estimates$pred, series=i))}
     #tframe(z) <- tframe(outputData(model))
     tfOnePlot(z,start=start,end=end, xlab=xlab[m+i], ylab=names[m+i]) # tsplot
     if(!is.null(Title) && (i==1) && (is.null(options()$PlotTitles)
                || options()$PlotTitles)) title(main = Title)
    }
  invisible()
  }
    


testEqual.TSestModel <- function(obj1, obj2, fuzz=0) # this could be better
  { testEqual.TSmodel( obj1$model, obj2$model, fuzz=fuzz) &
        testEqual.TSdata(obj1$data, obj2$data, fuzz=fuzz)
  }

testEqual.TSmodel <- function(obj1, obj2, fuzz=0)
{# return T if models are identical (excluding description)
  r       <- all(class(obj1) == class(obj2))
  if (r) r <-length(coef(obj1)) == length(coef(obj2))
  if (r) r <-all(fuzz >= abs(coef(obj1)   -  coef(obj2)))
  if (r) r <-length(obj1$location) == length(obj2$location)
  if (r) r <-all(obj1$location  ==     obj2$location)
  if (r) r <-length(obj1$i) == length(obj2$i)
  if (r) r <-all(obj1$i  ==     obj2$i)
  if (r) r <-length(obj1$j) == length(obj2$j)
  if (r) r <-all(obj1$j  ==     obj2$j)
  if (r) r <-length(obj1$const) == length(obj2$const)
  if (r) r <-all(obj1$const  ==     obj2$const)
  if (r) r <-length(obj1$const.location) == length(obj2$const.location)
  if (r) r <-all(obj1$const.location  ==     obj2$const.location)
  if (r) r <-length(obj1$const.i) == length(obj2$const.i)
  if (r) r <-all(obj1$const.i  ==     obj2$const.i)
  if (r) r <-length(obj1$const.j) == length(obj2$const.j)
  if (r) r <-all(obj1$const.j  ==     obj2$const.j)
  if (r)
    {if (is.ARMA(obj1))    r <-testEqual.ARMA(obj1,obj2, fuzz=fuzz)
     else if (is.SS(obj1)) r <-testEqual.SS(obj1,obj2, fuzz=fuzz)
    }
  r
}

testEqual.ARMA <- function(obj1, obj2, fuzz=0)
{    r <-length(obj1$l) == length(obj2$l)
     if (r) r <-all(obj1$l  ==     obj2$l)
     if (r) r <-length(obj1$const.l) == length(obj2$const.l)
     if (r) r <-all(obj1$const.l  == obj2$const.l)
     if (r) r <-length(obj1$A) == length(obj2$A)
     if (r) r <-all(fuzz >= abs(obj1$A   -     obj2$A))
     if (r) r <-length(obj1$B) == length(obj2$B)
     if (r) r <-all(fuzz >= abs(obj1$B   -     obj2$B))
     if (r) 
           {if (is.null(obj1$C)) r <- is.null(obj2$C)
             else
               {if (r) r <-length(obj1$C) == length(obj2$C)
                if (r) r <-all(fuzz >= abs(obj1$C   -     obj2$C))
           }  }
  r
}

testEqual.SS <- function(obj1,obj2, fuzz=0)
{    r <-length(obj1$F) == length(obj2$F)
     if (r) r <-all(fuzz >= abs(obj1$F   -     obj2$F))
     if (r) 
         {if(is.null(obj1$G)) r <- is.null(obj2$G)
          else
            {if (r) r <-length(obj1$G) == length(obj2$G)
             if (r) r <-all(fuzz >= abs(obj1$G   -     obj2$G))
         }  }
     if (r) r <-length(obj1$H) == length(obj2$H)
     if (r) r <-all(fuzz >= abs(obj1$H   -     obj2$H))
     if (is.innov.SS(obj1))
       {if (r) r <-length(obj1$K) == length(obj2$K)
        if (r) r <-all(fuzz >= abs(obj1$K   -     obj2$K))
       }
     else
       {if (r) r <-length(obj1$Q) == length(obj2$Q)
        if (r & (0 != length(obj2$Q)) ) r <-all(fuzz >= abs(obj1$Q - obj2$Q))
        if (r) r <-length(obj1$R) == length(obj2$R)
        if (r & (0 != length(obj2$R)) ) r <-all(fuzz >= abs(obj1$R - obj2$R))
       }
     if (r) if(is.null(obj1$z0)) r <- is.null(obj2$z0)
     else
       {if (r) r <-length(obj1$z0) == length(obj2$z0)
        if (r) r <-all(fuzz >= abs(obj1$z0   -     obj2$z0))
       }
     if (r) if(is.null(obj1$P0)) r <- is.null(obj2$P0)
     else
       {if (r) r <-length(obj1$P0) == length(obj2$P0)
        if (r) r <-all(fuzz >= abs(obj1$P0   -     obj2$P0))
       }
     if (r) if(is.null(obj1$rootP0)) r <- is.null(obj2$rootP0)
     else
       {if (r) r <-length(obj1$rootP0) == length(obj2$rootP0)
        if (r) r <-all(fuzz >= abs(obj1$rootP0   -     obj2$rootP0))
       }
  r
}


McMillanDegree <- function(model,  ...) UseMethod("McMillanDegree")

McMillanDegree.TSestModel <- function(model, ...)
 {McMillanDegree(TSmodel(model),...) }

McMillanDegree.ARMA <- function(model, fuzz=1e-4, verbose=TRUE, warn=TRUE, ...)
 {#  (... further arguments, currently disregarded)
    z  <- roots(model, warn=warn)
    gross <- length(z)
    zz <- outer(z,z,FUN="-")
    distinct <-sum(!apply((outer(1:length(z),1:length(z),FUN="<") & (Mod(zz) <fuzz)),2,any))
    deg <- list(gross=gross,distinct=distinct)
    if (verbose)
      {cat("Assuming the model is left prime:\n")
       cat("   Without distinguishing distinct roots the degree det(A(L)) =",deg$gross,"\n")
       cat("   Distinguishing  distinct  roots       the degree det(A(L)) =",deg$distinct,"\n")
       if(!is.null(model$TREND))
          {cat("The trend adds unit roots which are added to the degree. Multiple")
           cat("unit roots are not considered distinct (but probably should be).\n")
          }
      }
    invisible(deg)
 }

McMillanDegree.SS <- function(model, fuzz=1e-4, ...)
 {#  (... further arguments, currently disregarded)
   cat("state dimension = ", nstates(model),"/n")
   invisible()
 }

stability <- function(obj, fuzz=1e-4, eps=1e-15, digits=8, verbose=TRUE) UseMethod("stability")

stability.TSestModel <- function(obj, fuzz=1e-4, eps=1e-15, digits=8, verbose=TRUE)
   {stability(TSmodel(obj), fuzz=fuzz, eps=eps, digits=digits, verbose=verbose)}

stability.TSmodel <- function(obj, fuzz=1e-4, eps=1e-15, digits=8, verbose=TRUE)
  {stability(roots(obj, fuzz=fuzz, randomize=FALSE),
             fuzz=fuzz, eps=eps, digits=digits, verbose=verbose)}

stability.roots <- function(obj, fuzz=1e-4, eps=1e-15, digits=8, verbose=TRUE) 
   {#obj <- roots(model, fuzz=fuzz, randomize=FALSE)
    m <- Mod(obj)
    s <- if (all(m < (1.0 - eps))) TRUE else  FALSE
    if (verbose)
      {if (s) cat("The system is stable.\n")
       else   cat("The system is NOT stable.\n")
      }
    r <- cbind(obj, m)
    dimnames(r) <- list(NULL, c("Eigenvalues of F","moduli"))
    attr(s, "roots") <- r
    s
   }

stability.ARMA <- function(obj, fuzz=1e-4, eps=1e-15, digits=8, verbose=TRUE) 
   {z <- roots(obj, fuzz=fuzz, randomize=FALSE)
    if (is.null(z))
      {warning("The model has only a zero root.")
       return(TRUE)}
    m <- Mod(z)
    s <- if (all(m < (1.0 - eps))) TRUE else  FALSE
    if (verbose) {
        cat("Distinct roots of det(A(L)) and moduli are:\n")
        print(cbind(1/z, Mod(1/z)), digits = digits)
        if (!is.null(obj$TREND)) 
            cat("Trend not taken into account: ")
        if (s) 
            cat("The system is stable.\n")
        else cat("The system is NOT stable.\n")
      }
    r <- cbind(z, m)
    dimnames(r) <- list(NULL, c("Inverse of distinct roots of det(A(L))","moduli"))
    attr(s, "roots") <- r
    s
   }


roots <- function(obj, ...) UseMethod("roots")

roots.TSestModel <- function(obj, ...){roots(TSmodel(obj), ...)}

roots.SS <- function(obj, fuzz=0, randomize=FALSE, ...) 
{#  (... further arguments, currently disregarded)
    z <- eigen(obj$F, symmetric=FALSE, only.values=TRUE)$values
    if (randomize) if (sample(c(TRUE,FALSE),1)) z <- Conj(z)
      #this prevents + - ordering of complex roots (for Monte Carlo evaluations)
    classed(z,"roots")  # constructor (roots.SS)
}


roots.ARMA <- function(obj, fuzz=0, randomize=FALSE, warn=TRUE, by.poly=FALSE, ...) 
{#  (... further arguments, currently disregarded)
    if(1 == dim(obj$A)[1]) {
        warning("An ARMA model with no lags in the AR term has no roots.")
	return(NULL)
	}
    if(by.poly) z <- 1/polyrootDet(obj$A)
    else        z <- roots(toSS(obj))
    if (fuzz!=0)
      {zz <- outer(z,z,FUN="-")  # find distinct roots within fuzz
       z <- z[ !apply((outer(1:length(z),1:length(z),FUN="<")
                 & (Mod(zz) <fuzz)),2,any)]
      }
    # add unit roots for TREND elements.
    if (!is.null(obj$TREND))
      {#z <- c(rep(1,sum(0!=obj$TREND)), z)this depends on nature of TREND
       if (warn)
         warning("Unit roots may need to be added for non-zero trend elements.")
      }
    if (randomize) if (sample(c(TRUE,FALSE),1)) z <- Conj(z)
      #this prevents + - ordering of complex roots (for Monte Carlo evaluations)
    classed(z,"roots")  # constructor (roots.ARMA)
}


plot.roots <- function(x, pch='*', fuzz=0, ...)
{#  (... further arguments, currently disregarded)
 if (is.TSmodel(x))    x <- roots(x, fuzz=fuzz)
 if (is.TSestModel(x)) x <- roots(x, fuzz=fuzz)
 i <- 2*pi*(1:1000)/1000
 if (max(Mod(x)) <1 )
   {plot(sin(i),cos(i),pch='.', xlab='Im', ylab='Re')
    points(Re(x),Im(x), pch=pch)
   }
 else
   {#plot.default(Re(x),Im(x), pch=pch, xlab='Im', ylab='Re')
    plot(Re(x),Im(x), pch=pch, xlab='Im', ylab='Re')
    points(sin(i),cos(i),pch='.')
   }
 points(0,0,pch='+')
 invisible(x)
}


addPlotRoots <- function(v, pch='*', fuzz=0)
{if (is.TSmodel(v))    v <- roots(v, fuzz=fuzz)
 if (is.TSestModel(v)) v <- roots(v, fuzz=fuzz)
 points(Re(v),Im(v), pch=pch)
 invisible(v)
}


observability <- function(model)  UseMethod("observability")

observability.TSestModel <- function(model){observability(TSmodel(model)) }

observability.SS <- function(model){ 
  FF<-    model$F
  O <-    model$H
  HFn <- O
  for (n in 1:dim(FF)[1])  {
    HFn <- HFn %*% FF
    O <- rbind(O,HFn)
    }
  La.svd(O)$d
  }

observability.ARMA <- function(model){
  cat("not applicable to ARMA models\n")
  invisible()
  }


reachability <- function(model) UseMethod("reachability") 

reachability.TSestModel <- function(model){reachability(TSmodel(model))}

reachability.SS <- function(model){
 FF <-    model$F
 C  <-    model$G
 if (!is.null(C))
   {FnG <- C
    for (n in 1:dim(FF)[1])  
      {FnG <- FF %*% FnG
       C <- cbind(C,FnG)
      }
    cat("Singular values of reachability matrix for input: ",La.svd(C)$d)
   }
 if (is.innov.SS(model)) C <- model$K
 else      
  {C <- model$R
   if(dim(C)[1]==1)
     {if (any(C==0))
         {cat("State noise matrix is singular. All states are NOT excited!\n")
          return(C)
         }
      C <- 1/C
     }
   else
     {v<-La.svd(C)
      if (any(v$d==0))
         {cat("State noise matrix is singular. All states are NOT excited!\n")
          return(v$d)
         }
  #   C <-v$v %*% diag(1/v$d) %*% t(v$u) following is equivalent
  #   C <-v$v %*% (t(v$u) * 1/v$d) with svd (vs La.svd)
      C <-(v$u %*% v$vt ) * 1/v$d 
     }
   C <- model$Q %*% C
  }
 FnK <- C
 for (n in 1:dim(FF)[1])  
   {FnK <- FF %*% FnK
    C <- cbind(C,FnK)
   }
 d <- La.svd(C)$d
 cat("Singular values of reachability matrix for noise: ",d,"\n")
 invisible(d)
 }

reachability.ARMA <- function(model){ 
  cat("not applicable to ARMA models\n")
  invisible()
  }


checkBalance <- function(model) UseMethod("checkBalance") 

checkBalance.TSestModel <- function(model){checkBalance(TSmodel(model))}

checkBalance.SS <- function(model){ 
  FF<-    model$F
  O <-    model$H
  HFn <- O
  for (n in 1:dim(FF)[1])  {
    HFn <- HFn %*% FF
    O <- rbind(O,HFn)
    }
  #O <- t(O) %*% O 
  O <- crossprod(O) # observability gramian
  C <-    cbind(model$G,model$K)
  FnG <- C
  for (n in 1:dim(FF)[1])  {
    FnG <- FF %*% FnG
    C <- cbind(C,FnG)
    }
  C <- C %*% t(C) # controllability gramian
  difference <- O-C
  #cat("observability gramian minus controllability gramian:\n")
  #print(difference)
  cat("maximum absolute difference (O-C): ", max(abs(difference)),"\n")
  cat("maximum off-diagonal element of C: ", max(abs(C-diag(diag(C)))),"\n")
  cat("maximum off-diagonal element of O: ", max(abs(O-diag(diag(O)))),"\n")
  invisible()
  }

checkBalance.ARMA <- function(model){ 
  cat("not applicable to ARMA models\n")
  invisible()
  }


checkBalanceMittnik <- function(model) UseMethod("checkBalanceMittnik")

checkBalanceMittnik.TSestModel <- function(model)
   checkBalanceMittnik(TSmodel(model))

checkBalanceMittnik.SS <- function(model){ 
  FF  <-    model$F - model$K %*% model$H
  O   <-    model$H
  HFn <- O
  for (n in 1:dim(FF)[1])  {
    HFn <- HFn %*% FF
    O <- rbind(O,HFn)
    }
  #O <- t(O) %*% O 
  O <- crossprod(O)# observability gramian
  C <-    cbind(model$G,model$K)
  FnG <- C
  for (n in 1:dim(FF)[1])  {
    FnG <- FF %*% FnG
    C <- cbind(C,FnG)
    }
  C <- C %*% t(C) # controllability gramian
  difference <- O-C
  #cat("observability gramian minus controllability gramian:\n")
  #print(difference)
  cat("maximum absolute difference (O-C): ", max(abs(difference)),"\n")
  cat("maximum off-diagonal element of C: ", max(abs(C-diag(diag(C)))),"\n")
  cat("maximum off-diagonal element of O: ", max(abs(O-diag(diag(O)))),"\n")
  invisible()
  }

checkBalanceMittnik.ARMA <- function(model){ 
  cat("not applicable to ARMA models\n")
  invisible()
  }


############################################################

#     functions for model conversion   <<<<<<<<<<

############################################################

toSS <- function(model, ...) UseMethod("toSS")

toSS.TSestModel <- function(model, ...) 
	l(toSS(TSmodel(model), ...),TSdata(model))

toSS.SS <- function(model, ...) {model}
 #  (... further arguments, currently disregarded)
 
toSS.ARMA <- function(model,...){
    # convert an ARMA (or VAR) to a SS (innovations) representation
    if (is.null(model$A)) a<-0
    else a <- dim(model$A)[1] - 1  #order of polynomial arrays
    if (is.null(model$B)) b<-0
    else b <- dim(model$B)[1] - 1
    if (is.null(model$C)) cc<-0
    else cc<- dim(model$C)[1] - 1
    if ((b<=a) & (cc<=(a-1))) model <- toSSaugment(model)
    else  model <-toSSnested(model,...) #  (otherwise best working method) 
                  # A better approach would be an algorithm like Guidorzi's. 
    model
    }


toSSnested <- function(model, ...) UseMethod("toSSnested")

toSSnested.TSestModel <- function(model, ...) toSSnested(TSmodel(model), ...)

toSSnested.SS <- function(model, n=NULL, Aoki=FALSE, ...){
  #  (... further arguments, currently disregarded)
  # convert to a nested-balanced state space model by svd  a la Mittnik (or Aoki)
  if (is.null(n)) n <-ncol(model$F)  
  if (Aoki) stop("Aoki balancing no longer supported.") #return(Aoki.balance(model, n=n))
  else      return(balanceMittnik(model, n=n)) 
  }

toSSnested.ARMA <- function(model, n=NULL, Aoki=FALSE, ...){
  #  (... further arguments, currently disregarded)
  # convert to a nested-balanced state space model by svd  a la Mittnik (or Aoki)
  if (is.null(n)) n <- McMillanDegree(model)$distinct
  if (Aoki) stop("Aoki balancing no longer supported.") #return(Aoki.balance(model, n=n))
  else      return(balanceMittnik(model, n=n)) 
  }


toSSaugment <- function(model, ...) UseMethod("toSSaugment")

toSSaugment.TSestModel <- function(model, ...)
   l(toSSaugment(TSmodel(model), ...), TSdata(model))


toSSaugment.ARMA <- function(model, fuzz=1e-14, ...) {
  #  (... further arguments, currently disregarded)
  # convert by augmentation - state dimension may not be minimal
  # First sets A[1,,] = B[1,,] = I if that is not already the case.
   A <- model$A
   B <- model$B
   C <- model$C 
   if (fuzz  < max(abs(A[1,,]-diag(1,dim(A)[2]) )) )
      {A0.inv <- solve(A[1,,])
       A <-  polyprod(A0.inv,A)
       B <-  polyprod(A0.inv, B)
       if (!is.null(C)) C <- polyprod(A0.inv, C)
#       if (!is.null(TREND)) TREND <- t(A0.inv %*% t(TREND))
       }
   if (fuzz  < max(abs(B[1,,]-diag(1,dim(B)[2]) )) )
          B<- polyprod(solve(B[1,,]), B)
   if (is.null(model$A)) a<-0
   else a <- dim(model$A)[1] - 1  #order of polynomial arrays
   if (is.null(model$B)) b<-0
   else b <- dim(model$B)[1] - 1
   if (is.null(model$C)) cc<-0
   else cc<- dim(model$C)[1] - 1
   p <- nseriesOutput(model)          # Dim of endoenous Variables.
   m <-  nseriesInput(model)          # Dim of exogenous Variables.
   if (b>a) stop("The MA order cannot exceed the AR order to convert with state augmentation.")
   if (cc>(a-1)) stop(
      " The order of the input polynomial cannot exceed the AR order -1 to convert with state augmentation.")   
  #make three parameters A,B and C have convenient order by adding 0's.
   k <- 1 + a  
#  if (b != 0)
    {BB <- array(0,c(k,dim(B)[2:3]))
     BB[1:(b+1),,] <- B
    }
   if (m!=0) 
     {CC <- array(0,c(k,dim(C)[2:3]))  
      CC[1:(cc+1),,] <- C
     }
   FF <- matrix(NA,a*p,p)
   for (i in 1:a) FF[(1+p*(i-1)):(p*i),] <- -A[a-i+2,,]
   if(a>1) FF<-cbind(rbind(matrix(0,p,(a-1)*p),diag(1,(a-1)*p)),FF)
   if (m == 0) G <-NULL
   else
     {G <- matrix(NA,a*p,m)
      for (i in 1:a) G[(1+p*(i-1)):(p*i),] <- CC[a-i+1,,] 
     }
   H <- diag(1,p)
   if(a>1) H <- cbind( matrix(0,p,(p*(a-1))),H)
   K <- matrix(NA,a*p,p)
   for (i in 1:a) K[(1+p*(i-1)):(p*i),] <- -A[a-i+2,,]+BB[a-i+2,,]
   z0 <-NULL
   if(!is.null(model$TREND))   #add a constant state which feeds into the states
     {FF<-rbind(cbind(FF,0),0) # identified with outputs (through H).
      n <-dim(FF)[1]
      FF[n,n] <-1 
      if (p != length(model$TREND)) stop("This fails for matrix TREND.")
      FF[n-p:1,n] <- model$TREND
      z0 <- rep(0,n)
      z0[n] <-1
      H<-cbind(H,0)
      if (m!=0) G <- rbind(G,0)
      K<- rbind(K,0)
     }                     
   descr<-c(model$description,
            " Converted to state space by state augmentation.")
   SS(F.=FF,G=G,H=H,K=K,z0=z0,description=descr,
         input.names= seriesNamesInput(model),
        output.names=seriesNamesOutput(model))        
 }


gmap <- function(g, model) 
{# convert to an equivalent representation using a given matrix
 if (is.TSestModel(model)) model <- TSmodel(model)
 if  (!is.TSmodel(model)) stop("gmap expecting a TSmodel.")
 if ( is.SS(model))# transform State space model by g in GL(n)
  {n <- dim(model$F)[1]
   if (!is.matrix(g)) g<-diag(g,n) # if g is not a matrix make it into a diagonal matrix.
   if ((n !=dim(g)[1]) | (n !=dim(g)[2]) )
      stop("g must be a square matrix of dimensions equal the model state (or a scalar).")
   inv.g <- solve(g)
   model$F <-inv.g%*%model$F%*% g
   if (!is.null(model$G)) model$G <-inv.g %*%model$G
   model$H <-model$H %*% g
   if (!is.null(model$z0)) model$z0 <-c(inv.g %*%model$z0)
   if (is.innov.SS(model)) model$K <-inv.g %*% model$K
   if (is.nonInnov.SS(model)) 
      {model$Q <-inv.g %*% model$Q
       model$R <-model$R
      }       
 }
 if ( is.ARMA(model))
       {# if g is not a matrix make it into a diagonal matrix.
        if (! is.matrix(g)) g<- diag(g,dim(model$A)[2]) 
	for(l in 1:dim(model$A)[1]) model$A[l,  ,  ] <- g %*% model$A[l, ,]	
	for(l in 1:dim(model$B)[1]) model$B[l,  ,  ] <- g %*% model$B[l, ,]
	for(l in 1:dim(model$C)[1]) model$C[l,  ,  ] <- g %*% model$C[l, ,]
	if(!is.null(model$TREND))  model$TREND <- t(g %*% t(model$TREND))
       }
 setTSmodelParameters(model)
}


#findg <- function(model1,model2, minf=nlmin){ 
#  if (is.TSestModel(model1)[1]) model1 <- TSmodel(model1)
#  if   (!is.TSmodel(model1)) stop("findg expecting a TSmodel.")
#  if (is.TSestModel(model2)[1]) model2 <- TSmodel(model2)
#  if  (!is.TSmodel(model2 )) stop("findg expecting a TSmodel.")
#
#  if ( !is.SS(model1)| !is.SS(model2)) 
#      stop("findg only works for state space models")
#   n <- dim(model1$F)[1]
#   if ((n!= dim(model2$F)[1])
#     |(dim(model1$G)[2] != dim(model2$G)[2])
#     |(dim(model1$H)[1] != dim(model2$H)[1]))
#      stop("models must have the same dimensions for findg.")
#   para<- c(diag(1,n))
#   zzz.model1 <<- model1	  # This could be done with assign(frame=1  ??)
#   zzz.model2 <<- model2
#   zzz.n <<-n
#   func <- function(para){
#      gmodel1<- gmap(matrix(para,zzz.n,zzz.n),zzz.model1)
#      error <- 	gmodel1$F-zzz.model2$F
#      error <-c(error,(gmodel1$G-zzz.model2$G))
#      error <-c(error,(gmodel1$H-zzz.model2$H))
#      error <-c(error,(gmodel1$K-zzz.model2$K))
#      error <-c(error,(gmodel1$R-zzz.model2$R))
#      sum(error^2)
#   }
#   para <-minf(func,para)
#   rm(zzz.model1,zzz.model2,zzz.n)
#   matrix(para[[1]],n,n)
#   }


fixConstants <- function(model, fuzz=1e-5, constants=NULL){ 
  if (is.TSestModel(model)) model <- TSmodel(model)
  if  (!is.TSmodel(model)) stop("fixConstants expecting a TSmodel.")
  if (is.null(constants))
    {p <-abs(coef(model) - 1.0) < fuzz
     model$const <- c(model$const,rep(1.0,sum(p)))
     model$const.location <- c(model$const.location,model$location[p])
     model$const.i <- c(model$const.i,model$i[p])
     model$const.j <- c(model$const.j,model$j[p])
     if(is.ARMA(model)) model$const.l <- c(model$const.l,model$l[p])
     p <- (!p) & (abs(coef(model)) > fuzz) 
     model$coefficients <- coef(model)[p]
     model$location <- model$location[p]
     model$i <- model$i[p]
     model$j <- model$j[p]
     if(is.ARMA(model)) model$l <- model$l[p]
     return(setArrays(model))
    }
  else return(setTSmodelParameters(model,constants=constants))
  }


toSSinnov <- function(model, ...){
 #  (... further arguments, currently disregarded)
 data <- NULL
 if (is.TSestModel(model)) 
    {data <- TSdata(model)
     model <- TSmodel(model)
    }
 if  (!is.TSmodel(model)) stop("toSSinnov expecting a TSmodel.")
 if (!is.SS(model))  model <- toSS(model)
 if ( is.nonInnov.SS(model)) 
   {warning("this is not exactly correct.")
    PH  <-  model$Q %*% t(model$Q) %*% t(model$H)# use QQ' as state tracking error ??
    ft    <- (model$H %*% PH) + (model$R %*% t(model$R))        
    ft    <-  (ft + t(ft))/2   # force ft to be symmetric 
    model$K <-  t(solve(ft,t(model$F %*% PH))) 
    model$R <- NULL
    model$Q <- NULL
   }  
 model <- setTSmodelParameters(classed(model, c("innov","SS","TSmodel"))) # bypass constructor
 if (is.null(data)) model else l(model, data)
 }

toSSOform <- function(model) UseMethod("toSSOform")

toSSOform.TSestModel <- function(model) 
   l(toSSOform(TSmodel(model)), TSdata(model))

toSSOform.TSmodel <- function(model){
 if (!is.SS(model))       model <- toSS(model)
 if (!is.innov.SS(model)) model <- toSSinnov(model)
 n <- dim(model$H)[2]
 p <- nseriesOutput(model)
 if (p >= n) 
   {ginv <-model$H[1:n,]
    g <-solve(ginv)
   }
 else
   {sv   <- La.svd(model$H)
# svd    g    <- sv$v %*% diag(1/sv$d, ncol = length(sv$d))  %*%  t(sv$u) #right inv
# svd    ginv <- sv$u %*% diag( sv$d,  ncol = length(sv$d))  %*%  t(sv$v) #left  inv
    g    <- t(sv$vt) %*% diag(1/sv$d, ncol = length(sv$d))  %*%  t(sv$u) #right inv
    ginv <-   sv$u   %*% diag( sv$d,  ncol = length(sv$d))  %*%    sv$vt #left  inv
    g    <- cbind(g,   matrix(0,n, n-p))  # no good. these need to be full rank
    ginv <- rbind(ginv,matrix(0,n-p, n))  # and still convert H to [ I | 0 ]
    stop("This procedure is not working properly yet.") 
    # have fixed only nxp not yet nxn elements
   }
 model$F <- ginv %*% model$F %*% g
 if(!is.null(model$G)) model$G <- ginv %*% model$G
 model$H <- model$H %*% g
 model$K <- ginv %*% model$K  
 fixConstants(setTSmodelParameters(model))
}


fixF <- function(model){
  if (is.TSestModel(model)) model <- TSmodel(model)
  if  (!is.TSmodel(model)) stop("fixF expecting a TSmodel.")
  if (!is.SS(model))         model <- toSS(model)
  if ( is.nonInnov.SS(model))  model <- toSSinnov(model)
  p <-model$location == "f"
  model$const <- c(model$const, coef(model)[p])
  model$const.location <- c(model$const.location,model$location[p])
  model$const.i <- c(model$const.i,model$i[p])
  model$const.j <- c(model$const.j,model$j[p])
  p <- !p
  model$coefficients <- coef(model)[p]
  model$location <- model$location[p]
  model$i <- model$i[p]
  model$j <- model$j[p]
  cat("Remaining parameters: ",sum(p),"\n")
  if (is.null(model$G)) m <-0
  else  m <- dim(model$G)[2]
  n <- dim(model$F)[1]
  p <- dim(model$H)[1]
  cat("Theoretical parameter space dimension: ",n*(m+2*p),"\n")
  setArrays(model)
}


toSSChol <- function(model, ...) UseMethod("toSSChol")

toSSChol.TSestModel <- function(model, Om=NULL, ...) {
   #  (... further arguments, currently disregarded)
   if(is.null(Om)) Om <-model$estimates$cov
   l(toSSChol(TSmodel(model), Om=Om), TSdata(model))
   }

toSSChol.TSmodel <- function(model, Om=diag(1,nseriesOutput(model)), ...){
 #  (... further arguments, currently disregarded)
 if (!is.SS(model))  model <- toSS(model)
 if (is.innov.SS(model)) 
   {model$R <-t(chol(Om) )  # Om = RR'
    model$Q <- model$K %*% model$R
    model$K <- NULL
   }  
 classed(model, c( "nonInnov","SS","TSmodel" ) )  # bypass constructor
 }

toARMA <- function(model, ...) UseMethod("toARMA")

toARMA.TSestModel <- function(model, ...) 
	l(toARMA(TSmodel(model), ...), TSdata(model))

toARMA.ARMA <- function(model, ...) model

toARMA.SS <- function(model, fuzz=1e-10, ...){
    #  (... further arguments, currently disregarded)
    if (is.nonInnov.SS(model)) model <- toSSinnov(model)
    FF<-model$F
    G <-model$G
    H <-model$H
    K <-model$K
    m <-dim(G)[2]
    if (is.null(m)) m <-0
    n <-dim(FF)[1]
    p <-dim(H)[1]
    poly <-  - characteristicPoly(FF) # N.B. sign change  in Aoki vs Kailath
    if ( n != length(poly))
      stop("There is some problem. The characteristic polynomial length should = state dimension.")
    A <- array(NA,c(1+n,p,p))
    A[1,,] <-diag(1,p)
    for (i in 1:n) A[i+1,,] <- diag(-poly[i],p)
    if(any(is.na(A))) stop("error in calculation of A in toARMA.")

#                                       i
#            Fn [i,,] corresponds to   F    
    if (n > 1)
      {Fn <- array(NA,c(n-1,n,n))
       Fn[1,,] <- FF
      }
    if (n > 2) for (i in 2:(n-1))   Fn[i,,] <- FF %*% Fn[i-1,,]

    
    HFnK <- array(NA, c(n+1, p, p))
    HFnK[1, , ] <- diag(1, p)
    HFnK[2, , ] <- H %*% K
    if (n > 1) for (i in 3:(n+1)) HFnK[i, , ] <- H %*% Fn[i-2, , ] %*% K
    B <- array(NA, c(1+n, p, p))
    B[1, , ] <- diag(1, p)
    for (i in 1:n) 
       {B[i+1, , ] <- HFnK[i+1, , ]
        for (j in 1:i) B[i+1, , ] <- B[i+1, , ] - poly[j] *  HFnK[i+1-j, , ]
       }
    if(any(is.na(B))) stop("error in calculation of B in toARMA.")

    if (m == 0) C <- NULL
    else
      {C <- array(NA,c(n,p,m))   
       HFnG <- array(NA,c(n,p,m)) 
       HFnG[1,,] <- H %*% G
       if (n > 1) for (i in 2:n) HFnG[i,,] <- H %*% Fn[i-1,,] %*% G
       C[1,,] <- HFnG[1,,]
       for (i in 2:n)
         {C[i,,] <- HFnG[i,,]
          for (j in 1:(i-1)) C[i,,] <- C[i,,]-poly[j]*HFnG[i-j,,]
         }
       if(any(is.na(C))) stop("error in calculation of C in toARMA.")
      }
 fixConstants(ARMA(A=A,B=B,C=C, input.names =  seriesNamesInput(model),
                  output.names = seriesNamesOutput(model)), fuzz=fuzz)
}


#######################################################################

#                  Utility functions  <<<<<<<<<<

############################################################

#             functions for generating  statistics  
#        from data and for generating theoretical statistics
#           (i.e.- calculated from model not data)  

############################################################



Riccati <- function(A, B, fuzz=1e-10, iterative=FALSE)
{warning("This procedure has not been tested!")
 if (!iterative) 
  {n <- dim(A)[1]
   Atinv <-solve(t(A))    # symplectic matrix Vaughan (10)(12), R=0
   S   <- rbind(cbind(Atinv     , diag(0,n)),
               cbind(B %*% Atinv, A)) 
   Q <- eigen(S, symmetric=FALSE, only.values=FALSE)
   Q <- Q$vectors[,rev(order(Mod(Q$values)))]   # This may have imaginary parts.
   P <- Re( Q[(n+1):(n+n),1:n] %*% solve(Q[1:n,1:n]) ) #This should not have any significant im parts.
  }
 else
  {P<- diag(0,dim(A)[1])
   i <-0
   repeat    # yuk
     {P <- A%*% P %*% t(A) + B
      i <- i+1
      if (i>1000) break
      if (fuzz > max(abs(P-A %*% P %*% t(A) - B))) break
  }  }
 if (fuzz < max(abs(P-A %*% P %*% t(A) - B)))
      warning("Riccati failed! Result does not solve the Riccati equation!")
 P
}


markovParms <- function(model, blocks=NULL) 
{if (is.TSestModel(model)) model <- TSmodel(model)
 if  (!is.TSmodel(model)) stop("markovParms expecting a TSmodel.")
 if ( is.ARMA(model))
  #  M ={ Mi }={ Ci-1|Bi| -Ai}, i=2,...,k. k=max(a,b,cc). Assumes I=A[1,,]=B[1,,]
  {A <- model$A
   B <- model$B
   C <- model$C
   a <- dim(A)[1] - 1      # order of polynomial arrays
   if (is.na(a)) a <- 0
   b <- dim(B)[1] - 1
   if (is.na(b)) b <- 0
   cc <- dim(C)[1] - 1 
   if (is.na(cc)) cc <- 0
   m <- dim(C)[3]          # Dim of exogenous Variables.
   if (is.null(m))  m <- 0                         
   #make three parameters A,B and C have convenient order by adding 0's.    
   if (is.null(blocks))
      blocks <- 1 + max(2,a,cc,b) 
    # if blocks is not at least 3 the Hankel shift does not work
   if (a != 0) 
    {AA <- array(0,c(blocks,dim(A)[2:3]))
     AA[1:(a+1),,] <- A
    }
   if (b != 0)
    {BB <- array(0,c(blocks,dim(B)[2:3]))
     BB[1:(b+1),,] <- B
    }
   if (m != 0)
    {CC <- array(0,c(blocks,dim(C)[2:3]))  
     CC[1:(cc+1),,] <- C
    }
   M <- NULL
   if (b != 0) cat (" Warning: This has only been developed for SS and VARX models.\n")
   if(!is.null(model$TREND)) cat(" Warning: Markov parameter generation does not account for trends.")
   for(i in 2:blocks)  
       {if (m != 0) M <- cbind(M,CC[(i-1),,]) 
        if (b != 0) M <- cbind(M, BB[i,,])      # constant term ignored (=I)
        if (a != 0) M <- cbind(M,-AA[i,,])      # constant term ignored (=I)
       }
  }   
else if ( is.SS(model))
   {FF<-model$F
    G <-model$G
    H <-model$H
    if (is.innov.SS(model)) K <-model$K
    else               K <- model$Q %*% solve(model$R)
    if (is.null(blocks)) blocks <- 1+dim(FF)[1]
    FF <- FF - K %*% H  # model transformed a la Mittnik so lagged outputs are inputs
    FnGK <- cbind(G,K)
    M <- H %*% FnGK
    i<-0 # M should have at least 2 blocks or Hankel shift does not work.
    stopp <- FALSE
    while((i<=blocks) & !stopp)   #  no. of blocks affect Hankel size
       {FnGK <- FF %*% FnGK
        M <- cbind(M,H %*% FnGK)
        i <-i+1   # count should not be necessary, but insures an end.
        stopp <- (i>3) & ( max(abs(FnGK)) < 1e-15)
       }
  }
else stop("markovParms requires an ARMA or SS model.")
M       
}



############################################################

#     polynomial utility functions   <<<<<<<<<<

############################################################


polyprod <- function(a,b)
{ # product of two polynomials.
# The convention used is by polyvalue and polyroot is constant first,
#  highest order coef. last. The reverse convention could also be used for multiplication.
# This function handles scalar (ie. non-matrix) and matrix polynomials.
# Scalar polynomials are vectors of length 1+the polynomial order.
# Polynomial matrices are defined as 3 dimensional arrays with the last 2
# dimensions as the matrix dimension and the first equal 1+the
# polynomial order.

    pprod <- function(a,b)  # local function, product of non-matrix polys.
       {n <- length(a) +length(b) -1
        if (is.null(a))    return(NA)
        if (is.null(b))    return(NA)
        if (0 ==length(a)) return(NA)
        if (0 ==length(b)) return(NA)
        if (any(is.na(a))) return(NA)
        if (any(is.na(b))) return(NA)
	r <- rep(NA, n)
        z <- outer(a, b) 
        zi <- 1 + outer(1:length(a)-1,1:length(b)-1,"+")
	for(i in 1:n) r[i]<- sum(z[i==zi])
	r
       }
   psum <- function(a,b)  # local function, sum of non-matrix polys.
    {if (length(a) < length(b)) return(c(a,rep(0,length(b)-length(a))) + b)
     else                       return(c(b,rep(0,length(a)-length(b))) + a)
    }
   if (is.vector(b)  && (is.array(a) | is.matrix(a)))
     {z <- b; b <- a; a <- z } # scalar multiplication commutes (even for scalar polynomials)
   if (is.null(a))      r <- NULL  
   else if (is.null(b)) r <- NULL   
   else if (is.vector(a))  
          {if      (is.vector(b)) r <-pprod(a,b)
           else if (is.matrix(b))
              {r <- array(NA,c(length(a),dim(b)))
               for (i in 1:(dim(b)[1])) 
                  for (j in 1:(dim(b)[2]))
                     r[,i,j] <- pprod(a,b[i,j])
              }
           else if (is.array(b))
              {r <- array(NA,c(length(a)+dim(b)[1]-1,dim(b)[2:3]))
               for (i in 1:(dim(b)[2])) 
                  for (j in 1:(dim(b)[3]))
                     r[,i,j] <- pprod(b[,i,j],a)
              }
          }
   else if (is.matrix(a))
          {if (is.matrix(b)) r <- a %*% b
           else if (is.array(b))
              {if (dim(a)[2] != dim(b)[2]) 
                  stop("Matrix polynomial dimensions do not conform.")
               r <- array(0,c(dim(b)[1],dim(a)[1], dim(b)[3]))
               #for (i in 1:(dim(a)[1])) for (j in 1:(dim(b)[3]))
                 #for (k in 1:(dim(a)[2])) r[,i,j] <- psum(r[,i,j], pprod(a[i,k],b[,k,j]))
               for (k in 1:(dim(b)[1]))
                 r[k,,] <- a %*% array(b[k,,],dim(b)[2:3])
                 #array above insures b a matrix, drop=FALSE cannot be used
              }
          }
   else if (is.array(a))
        if (is.matrix(b))
          {if (dim(a)[3] != dim(b)[1])
              stop("Matrix polynomial dimensions do not conform.")
           r <- array(0,c(dim(a)[1],dim(a)[2], dim(b)[2]))
           #for (i in 1:(dim(a)[2])) for (j in 1:(dim(b)[2]))
           #  for (k in 1:(dim(a)[3]))  r[,i,j] <- psum(r[,i,j], pprod(a[,i,k],b[k,j]))
           for (k in 1:(dim(a)[1]))
                 r[k,,] <- array(a[k,,],dim(a)[2:3]) %*% b
                 #array above insures b a matrix, drop=FALSE cannot be used
          }
        else if (is.array(b))
          {if (dim(a)[3] != dim(b)[2]) 
              stop("Matrix polynomial dimensions do not conform.")
           r <- array(0,c(dim(b)[1]+dim(a)[1]-1,dim(a)[2], dim(b)[3]))
           for (i in 1:(dim(a)[2])) 
              for (j in 1:(dim(a)[3]))
                 for (k in 1:(dim(a)[3]))
                   r[,i,j] <- psum(r[,i,j], pprod(a[,i,k],b[,k,j]))
          }
   else stop("polynomial product not defined for these objects")
r
}



polysum <- function(a,b){
   # sum of two polynomials (including polynomial arrays)
    psum <- function(a,b)  # local function for non-matrix polys.
 {if (length(a) < length(b))  {r <- b;  r[1:length(a)] <-r[1:length(a)] + a }
  else  {r <- a; r[1:length(b)] <-r[1:length(b)] + b }
  r 
 }
   if (is.null(a) )     r <- b
   else if (is.null(b)) r <- NULL
   else if (is.vector(a) && is.vector(b)) r <-psum(a,b)
   else if (is.array(a) )
        if (is.array(b))
           if ( all(dim(a)[2:3] == dim(b)[2:3]))
             {if (dim(a)[1] < dim(b)[1])
                    {r <- b
                     r[1:(dim(a)[1]),,] <-r[1:(dim(a)[1]),,] + a
                    }
                 else
                    {r <-  a
                     r[1:(dim(b)[1]),,] <-r[1:(dim(b)[1]),,] + b
                    }
             }
           else stop("polynomial matrix dimensions must agree")
        else if (is.vector(b))
          {r <- array(NA,c(max(length(b),dim(a)[1]),dim(a)[2:3]))
           for (i in 1:(dim(a)[2])) 
              for (j in 1:(dim(a)[3]))
                 r[,i,j] <- psum(a[,i,j],b)
          }
   else if (is.vector(a) && is.array(b))
          {r <- array(NA,c(max(length(a),dim(b)[1]),dim(b)[2:3]))
           for (i in 1:(dim(b)[2])) 
              for (j in 1:(dim(b)[3]))
                 r[,i,j] <- psum(b[,i,j],a)
          }
   else stop("polynomial sum not defined for these objects")
r             
}

polyrootDet <- function(a)
{#roots of the determinant of a.  Note: polydet is slow. There is room for improvement here!
z <- polydet(a)
if (length(z)==1) 
  stop(paste("root cannot be calculated for degree 0 determinant = ", as.character(z)))
polyroot(z)
}

polydet <- function(a)
{# recursive pessimist: Life is just one damned thing before another.
 # attributed to Richard Bird in The Mathematical Intelligencer, Winter 1994.
 #Recursively form the determinant of a polynomial matrix a, where the first
 #  dim stores the coefs. of the polynomial for each a[,i,j].
	n <- dim(a)[2]  
        if (n != dim(a)[3])
           stop( "The determinant is only defined for square polynomial matrices")
	if(1 == n) r <- c(a)
	else 
          {r<- 0   # previously NULL
           for (i in 1:n) 
             {if(!all(0==a[,i,1]))#not nec.but faster for sparse arrays
                   r<- polysum(r,(-1)^(i+1)*
          polyprod(a[,i,1],Recall(a[,(1:n)[i!=(1:n)],2:n,drop=FALSE])))
              r[is.na(r)] <- 0
              if (any(0==r)) 
                   {if (all(r==0)) r <- 0
                    else  r <-r[0 == rev(cumprod(rev(r==0)))] #remove trailing zeros
                    if(0==length(r)) r <- 0
             }     } 
          }
r
}

polyvalue <- function(coef, z)
{# evaluate a polynomial given by coef (constant first) at z
 # could be extended for matrix coef.
#           n-1           n-2
#  coef[n]*z   + coef[n-1]*z   + ... + coef[1]  
  coef %*% z^(0:(length(coef)-1))
}


characteristicPoly <- function(a){
	# coefficients of the characteristic polynomial of a matrix
	# return a vector of the coefficients of the characteristic polynomial of a.
	# ref. Kailath "Linear Systems" p657

   tr <- function(a){ # calc the trace of a matrix
     sum(diag(a))
   }
   n <- dim(a)[1]
   if (n != dim(a)[2]) stop(" arguement must be a square matrix.")
   s <-array(0,c(n,n,n))
   if (n==1) a2 <- -a
   else
     {a2 <- rep(0,n)  
      s[1,,] <- diag(1,n)
      for (i in 1:(n-1))
        {a2[i] <- -tr(s[i,,]%*%a)/i
         s[i+1,,] <- (s[i,,]%*%a) + diag(a2[i],n)
        }
      a2[n] <- -tr(s[n,,]%*%a)/n
     }
   a2
}
companionMatrix <- function(a)
{# return the (top) companion matrix for a 3 dim array a (polynomial matrix), 
#  where the 1st dim corresponds to coefs. of polynomial powers (lags).
# ref. Kailath "Linear Systems" p659
  p <- dim(a)[2]
  if (p!= dim(a)[3]) 
    stop("companion matrix can only be computed for square matrix polynomials")
  l <- dim(a)[1]  # 1+ order of a
  C <- rbind(matrix(0,p,l*p), cbind(diag(1,(l-1)*p,(l-1)*p), matrix(0,(l-1)*p,p)))
  for (i in 1:l) C[1:p,((l-1)*p+1):(l*p)] <- -a[l,,]
  C
}

############################################################

#     internal utility functions    <<<<<<<<<<

############################################################

nstates <- function(x) UseMethod( "nstates")
nstates.SS <- function(x)  nrow(x$F) 
nstates.ARMA <- function(x) stop("ARMA models do not have a state vector.")
nstates.TSestModel <- function(x) nstates(TSmodel(x))

nseriesInput <- function(x) UseMethod( "nseriesInput")
nseriesInput.default <- function(x)
   if (is.null(x$input)) 0 else nseries(x$input)
nseriesInput.SS <- function(x)   {if (is.null(x$G)) 0 else  dim(x$G)[2] }
nseriesInput.ARMA <- function(x) {if (is.null(x$C)) 0 else dim(x$C)[3] }
nseriesInput.TSestModel <- function(x){nseriesInput(x$data)}

nseriesOutput <- function(x)UseMethod("nseriesOutput")
nseriesOutput.default <- function(x)
   if (is.null(x$output)) 0 else nseries(x$output)
nseriesOutput.SS <- function(x){dim(x$H)[1] }
nseriesOutput.ARMA <- function(x){dim(x$A)[2] }
nseriesOutput.TSestModel <- function(x){nseriesOutput(x$data)}


checkConsistentDimensions <- function(obj1, obj2=NULL)
   # check data & model dimensions
   UseMethod("checkConsistentDimensions")
   
checkConsistentDimensions.TSdata <- function(obj1, obj2=NULL)
  {if(is.null(obj2))
     stop("A TSmodel must be supplied as the second argument to this method.")
   checkConsistentDimensions(obj2, obj1) #(data,model) -> (model, data)
   }
checkConsistentDimensions.TSestModel <- function(obj1, obj2=NULL)
   {if(is.null(obj2)) obj2 <- obj1$data
    checkConsistentDimensions(obj1$model, obj2)
   }

checkConsistentDimensions.SS <- function(obj1, obj2=NULL) # (model, data)
 {m <- ncol(obj1$G)
  n <- nrow(obj1$F)
  p <- nrow(obj1$H)
  if (n!= ncol(obj1$F)) stop("Model F matrix must be square.")
  if (n!= ncol(obj1$H))
      stop("Model H matrix have second dimension consistent with matrix F.")
  if (!is.null(obj1$G)) if(n!= nrow(obj1$G))
      stop("Model G matrix have first dimension consistent with matrix F.")
  if (!is.null(obj1$K)) if(n!= nrow(obj1$K))
      stop("Model K matrix have first dimension consistent with matrix F.")
  if (!is.null(obj1$K)) if(p!= ncol(obj1$K))
      stop("Model K matrix have second dimension consistent with matrix H.")
  if (!is.null(obj1$Q)) if(n!= nrow(obj1$Q))
      stop("Model Q matrix must have first dimension consistent with matrix F.")
  if (!is.null(obj1$R)) if(p!= nrow(obj1$R))
      stop("Model R matrix must have first dimension consistent with matrix H.")
  if (!is.null(obj1$R)) if(p!= ncol(obj1$R))
      stop("Model R matrix must be square.")

  if (!is.null(obj2))
   {if(dim(obj1$H)[1] != nseriesOutput(obj2))
       stop("Model and data output dimensions do not correspond.\n")
    if(is.null(obj1$G))
      {if(0 != nseriesInput(obj2))
        stop("Model and data input dimensions do not correspond.\n")
      }
    else
      {if(dim(obj1$G)[2] != nseriesInput(obj2))
         stop("Model and data input dimensions do not correspond.\n")
      }
   }
  return(TRUE)
 }

checkConsistentDimensions.ARMA <- function(obj1, obj2=NULL) #(model, data)
 {p <-dim(obj1$A)[2]
  if (p!= dim(obj1$A)[3]) stop("Model A array dim 2 and 3 should be equal.")
  if (p!= dim(obj1$B)[2]) stop("Model B array dim inconsistent with array A.")
  if (p!= dim(obj1$B)[3]) stop("Model B array dim inconsistent with array A.")
  if (!is.null(obj1$C))
      if (p!= dim(obj1$C)[2]) stop("Model C array dim inconsistent with array A.")

  if (!is.null(obj2))
   {if(dim(obj1$A)[2] != nseriesOutput(obj2))
       stop("Model and data output dimensions do not correspond.\n")
    if(is.null(obj1$C))
      {if(0 != nseriesInput(obj2))
        stop("Model and data input dimensions do not correspond.\n")
      }
    else
      {if(dim(obj1$C)[3] != nseriesInput(obj2))
         stop("Model and data input dimensions do not correspond.\n")
      }
   }
  return(TRUE)
 }

checkConsistentDimensions.default <- function(obj1, obj2=NULL)
   {stop(paste("No method for obj1ect of class ", class(obj1), "\n"))}



TSestModel <- function(obj) UseMethod("TSestModel")
TSestModel.TSestModel <- function(obj) {obj} # return everything


TSmodel <- function(obj, ...) UseMethod("TSmodel")
TSmodel.TSmodel <- function(obj, ...){obj}
#  (... further arguments, currently disregarded)

TSmodel.TSestModel <- function(obj, ...)
  {#  (... further arguments, currently disregarded)
   # Return a TSmodel object but also retains convergence info (if not null).
   obj$model$converged <- obj$converged
   obj$model
  }




ARMA <- function(A=NULL, B=NULL, C=NULL, TREND=NULL, 
          constants=NULL,
	  description=NULL, names=NULL, input.names=NULL, output.names=NULL) 
  {if  (is.null(A)) stop("specified structure is not correct for ARMA model.")
   # and fix some simple potential problems
   if(is.null(dim(A)))   A <- array(A, c(length(A),1,1))
   if(is.null(B)) stop("B array must be specified for ARMA class models.")
   if(is.null(dim(B)))   B <- array(B, c(length(B),1,1))
   if(2==length(dim(B))) B <- array(B, c(1, dim(B)))
   if(!is.null(C) && is.null(dim(C))) C <- array(C, c(length(C),1,1))
   model <- list(A=A, B=B, C=C, TREND=TREND, 
                 constants=constants, description=description)
   class(model) <- c("ARMA","TSmodel") # constructor
   if(!is.null(names)) seriesNames(model) <- names
   else
     {if(!is.null( input.names))  seriesNamesInput(model) <- input.names
      if(!is.null(output.names)) seriesNamesOutput(model) <- output.names
     }
   checkConsistentDimensions(model)
   setTSmodelParameters(model)
  }



SS <- function(F.=NULL, G=NULL, H=NULL, K=NULL, Q=NULL, R=NULL,
          z0=NULL, P0=NULL, rootP0=NULL,
          constants=NULL,
          description=NULL, names=NULL, input.names=NULL, output.names=NULL)   
  {if (is.null(F.) | is.null(H))
       stop("specified stucture is not correct for SS model.")
   # and fix some simple potential problems
   if(1==length(F.))  
     {if(!is.matrix(F.))                 F. <- matrix(F.,1,1)
      if(!is.matrix(H))                  H  <- matrix(H,length(H),1)
      if(!is.null(G) && !is.matrix(G))   G  <- matrix(G,1,length(G))
      if(!is.null(K) && !is.matrix(K))   K  <- matrix(K,1,length(K))
      if(!is.null(Q) && !is.matrix(Q))   Q  <- matrix(Q,1,length(Q))
     }
   model <- list(F=F., G=G, H=H, K=K, Q=Q, R=R, z0=z0, P0=P0, rootP0=rootP0,
                 constants=constants, description=description)
   if (!is.null(model$K)) class(model) <- c("innov","SS","TSmodel" ) else
   if (!is.null(model$Q)) class(model) <- c( "nonInnov","SS","TSmodel") # constructor
   else stop("specified stucture is not correct for SS model.")
   if(!is.null(names)) seriesNames(model) <- names
   else
     {if(!is.null( input.names))  seriesNamesInput(model) <-  input.names
      if(!is.null(output.names)) seriesNamesOutput(model) <- output.names
     }
   checkConsistentDimensions(model)
   setTSmodelParameters(model)
  }



"coef<-" <- function(object, value) UseMethod("coef<-")
"coef<-.default" <- function(object, value) {
   object$coefficients <- value
   object
   }
   
coef.TSmodel <- function(object, ...) object$coefficients
#  (... further arguments, currently disregarded)

coef.TSestModel <- function(object, ...) coef(TSmodel(object))
#  (... further arguments, currently disregarded)

is.TSmodel <- function(obj){inherits(obj,"TSmodel")}
is.TSestModel <- function(obj){inherits(obj,"TSestModel")}
is.SS <- function(obj){inherits(obj,"SS")}
is.innov.SS <- function(obj){inherits(obj,"SS")& inherits(obj,"innov")}
is.nonInnov.SS <- function(obj){inherits(obj,"SS")&inherits(obj,"nonInnov")}
is.ARMA <- function(obj){inherits(obj,"ARMA")}

# complete parameter info. based on representation info. 

setTSmodelParameters <- function(model, constants=model$constants)  
   UseMethod("setTSmodelParameters")

setTSmodelParameters.TSestModel <- function(model, constants=TSmodel(model)$constants)  
  setTSmodelParameters(TSmodel(model), constants=constants)


setTSmodelParameters.SS <- function(model, constants=model$constants) { 

 locateSS <- function(A,Ac,a,I,J,plist)# local function for locating parameters
  {indicate <-  (A==1.0)                # constants
   if (!is.null(Ac)) indicate <- indicate | Ac  # Ac is T for fixed entries
   plist$const <- c(plist$const,A[indicate])
   plist$const.location <- c(plist$const.location,rep(a,sum(indicate)))
   if (is.vector(A))  # z is not a matrix, (a=="z") also should work
     {plist$const.i <- c(plist$const.i,(1:length(A))[indicate]) 
      plist$const.j <- c(plist$const.j,(1:length(A))[indicate]) #dup.but ok
     }
   else if(is.matrix(A))
     {plist$const.i <- c(plist$const.i,row(A)[indicate]) 
      plist$const.j <- c(plist$const.j,col(A)[indicate])
     }
   else stop("The dimension of something in the SS model structure is bad.")
   indicate <- (A!=0.0) & (A!=1.0) 
   if (!is.null(Ac)) indicate <- indicate & (!Ac)
   coef(plist) <- c(coef(plist), A[indicate])          # parameters
   plist$location <- c(plist$location,rep(a,sum(indicate)))
   if (is.vector(A))  # z is not a matrix, (a=="z") also should work
     {plist$i <- c(plist$i,(1:length(A))[indicate]) 
      plist$j <- c(plist$j,(1:length(A))[indicate]) #dup.but ok
     }
   else if(is.matrix(A))
     {plist$i <- c(plist$i,row(A)[indicate])
      plist$j <- c(plist$j,col(A)[indicate])
     }
   else stop("The dimension of something in the SS model structure is bad.")
   plist
 }# end locateSS    

    m <-dim(model$G)[2]
    n <-dim(model$F)[1]
    if (n!= dim(model$F)[2]) stop("Model F matrix must be square.")
    if (n!= dim(model$H)[2])
      stop("Model H matrix have second dimension consistent with matrix F.")
    if (!is.null(model$G)) if(n!= dim(model$G)[1])
      stop("Model G matrix have first dimension consistent with matrix F.")
    p <-dim(model$H)[1]
    plist <- list(coefficients=NULL,location=NULL,i=NULL,j=NULL,
                 const=NULL,const.location=NULL,const.i=NULL,const.j=NULL)
    plist <- locateSS(model$F, constants$F,"f",n,n,plist)
    if(!is.null(m)) plist <- locateSS(model$G,constants$G,"G",n,m,plist)
    plist <- locateSS(model$H,constants$H,"H",p,n,plist)
    if(!is.null(model$z0)) plist <- locateSS(model$z0,constants$z0,"z",p,n,plist)
    if(!is.null(model$P0)) plist <- locateSS(model$P0,constants$P0,"P",p,n,plist)
    if(!is.null(model$rootP0)) plist <- locateSS(model$rootP0,constants$rootP0,"rP",p,n,plist)
    if (is.innov.SS(model)) 
      {plist <- locateSS(model$K,constants$K,"K",n,p,plist)
       # note constants$H, etc are logical arrays (if not NULL) to indicate
       #      parameters which are to remain fixed, so that setTSmodelParameters
       #       knows to put them in const. This feature has not been used much.
      }
    else
      {plist <- locateSS(model$Q,constants$Q,"Q",n,n,plist)
       plist <- locateSS(model$R,constants$R,"R",p,p,plist)
      }

#    list.add(model, names(plist) ) <- plist
    model[ names(plist) ] <- plist
    model
}

setTSmodelParameters.ARMA <- function  (model, constants=model$constants) { 

 locateARMA <- function(A,Ac, a,I,J,L,plist){ # local function for locating parameters
  ind <- function(x, i) # equivalent of row and col for higher dim arrays.
   {
	d <- dim(x)
	id <- 1:length(d)
	id <- c(i, id[i!=id])
	y <- array(1:(d[i]), dim(x)[id])
	aperm(y, order(id))
   }
   indicate <-  (A==1.0)                # constants
   if (!is.null(Ac)) indicate <- indicate | Ac  # Ac is T for fixed entries
   plist$const <- c(plist$const,A[indicate])
   plist$const.location <- c(plist$const.location,rep(a,sum(indicate)))
   if (is.vector(A))# trend is a vector not an array (a=="t")
     {plist$const.l <- c(plist$const.l, rep(0,sum(indicate))) 
      plist$const.i <- c(plist$const.i, (1:length(A))[indicate]) 
      plist$const.j <- c(plist$const.j, rep(0,sum(indicate)))
     }
   else if (is.matrix(A))  # trend is a matrix
     {plist$const.l <- c(plist$const.l, rep(0,sum(indicate))) 
      plist$const.i <- c(plist$const.i, row(A)[indicate]) 
      plist$const.j <- c(plist$const.j, col(A)[indicate])
     }
   else if(3 == length(dim(A)))
     {plist$const.l <- c(plist$const.l, ind(A,1)[indicate]) 
      plist$const.i <- c(plist$const.i, ind(A,2)[indicate]) 
      plist$const.j <- c(plist$const.j, ind(A,3)[indicate])
     }
   else stop("The dimension of something in the ARMA model structure is bad.")
   indicate <- (A!=0.0) & (A!=1.0)
   if (!is.null(Ac)) indicate <- indicate & (!Ac)
   coef(plist) <- c(coef(plist), A[indicate])          # parameters
   plist$location <- c(plist$location,rep(a,sum(indicate)))
   if (is.vector(A))# trend is a vector not an array (a=="t")
     {plist$l <- c(plist$l, rep(0,sum(indicate))) 
      plist$i <- c(plist$i, (1:length(A))[indicate]) 
      plist$j <- c(plist$j, rep(0,sum(indicate)))
     }
   else if (is.matrix(A))  # trend is a matrix
     {plist$l <- c(plist$l, rep(0,sum(indicate))) 
      plist$i <- c(plist$i, row(A)[indicate]) 
      plist$j <- c(plist$j, col(A)[indicate])
     }
   else if(3 == length(dim(A)))
     {plist$l <- c(plist$l, ind(A,1)[indicate] )
      plist$i <- c(plist$i, ind(A,2)[indicate] )
      plist$j <- c(plist$j, ind(A,3)[indicate])
     }
   else stop("The dimension of something in the ARMA model structure is bad.")
   plist  
 }# end locateARMA

       m <-dim(model$C)[3]
       p <-dim(model$A)[2]
       a <-dim(model$A)[1]
       b <-dim(model$B)[1]
       cc <-dim(model$C)[1]
       if (p!= dim(model$A)[3]) stop("Model A array dim 2 and 3 should be equal.")
       if (p!= dim(model$B)[2]) stop("Model B array dim inconsistent with array A.")
       if (p!= dim(model$B)[3]) stop("Model B array dim inconsistent with array A.")
       if (!is.null(model$C))
          if (p!= dim(model$C)[2]) stop("Model C array dim inconsistent with array A.")

       plist <- list(coefficients=NULL,location=NULL,i=NULL,j=NULL,
                         const=NULL,const.location=NULL,
                         const.i=NULL,const.j=NULL,l=NULL,const.l=NULL)
       plist <- locateARMA(model$A, constants$A, "A",p,p,a,plist)
       plist <- locateARMA(model$B, constants$B, "B",p,p,b,plist)
       if(!is.null(cc)) plist <- 
                locateARMA(model$C, constants$C, "C",p,m,cc,plist)
       if(!is.null(model$TREND)) plist <- 
	        locateARMA(model$TREND, constants$TREND, "t",p,m,cc,plist)

#       list.add(model, names(plist) ) <- plist
       model[ names(plist) ] <- plist
       model
}


setArrays <- function(model, coefficients=NULL, constants=NULL)  
  # complete representaion info. based on parameter info.  
  UseMethod("setArrays")
   
setArrays.TSestModel <- function(model, coefficients=NULL, constants=NULL)  
 setArrays(TSmodel(model), coefficients=coefficients, constants=constants) 
    
setArrays.SS <- function(model, coefficients=NULL, constants=NULL){
	# N.B. Dimension and class (innov/ nonInnov) info. is assumed accurate
    if (is.null(coefficients)) coefficients   <- coef(model)
                        else   model$coefficients <- coefficients
    a.pos  <- model$location
    i.pos  <- model$i
    j.pos  <- model$j
    const  <-if (is.null(constants)) model$const else constants #untested
    ca.pos <- model$const.location
    ci.pos <- model$const.i
    cj.pos <- model$const.j  
    m <-dim(model$G)[2]
    n <-dim(model$F)[1]
    p <-dim(model$H)[1]
    f       <-  matrix(0,n,n)    # F     
    if(!is.null(m)) G        <-  matrix(0,n,m)    # control 
    H       <-  matrix(0,p,n)    # H       
    K        <-  matrix(0,n,p)    # Kalman gain
    if (is.null(model$Q)) Q <-  matrix(0,n,n)        # eventually discarded
    else                  Q <- array(0,dim(model$Q)) # system noise might not be nxn
    R        <-  diag(0,p,p)     # measurement noise
    z       <-  rep(0,n)        # initial state
    P       <-  diag(0,n)       # initial tracking error
    rP      <-  diag(0,n)       # root initial tracking error
    if (length(coefficients)>0) 
       {i <- a.pos == "f"
        f[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        if(!is.null(m)) 
          {i <- a.pos == "G"
           G[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
          }
        i <- a.pos == "H"
        H[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        i <- a.pos == "K"
        K[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        i <- a.pos == "Q"
        Q[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        i <- a.pos == "R"
        R[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        i <- a.pos == "z"
        z[i.pos[i]] <- coefficients[i]
        i <- a.pos == "P"
        P[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
        i <- a.pos == "rP"
        rP[cbind(i.pos[i],j.pos[i])] <- coefficients[i]
       }
    if (length(const)>0) 
       {i <- ca.pos == "f"
        f[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        if(!is.null(m)) 
          {i <- ca.pos == "G"
           G[cbind(ci.pos[i],cj.pos[i])] <- const[i]
          }
        i <- ca.pos == "H"
        H[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        i <- ca.pos == "K"
        K[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        i <- ca.pos == "Q"
        Q[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        i <- ca.pos == "R"
        R[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        i <- ca.pos == "z"
        z[ci.pos[i]] <- const[i]
        i <- ca.pos == "P"
        P[cbind(ci.pos[i],cj.pos[i])] <- const[i]
        i <- ca.pos == "rP"
        rP[cbind(ci.pos[i],cj.pos[i])] <- const[i]
       }
    model$F <- f 
    if(!is.null(m)) model$G <- G
    model$H <- H
    if (is.innov.SS(model))
      model$K <- K  
    else
      {model$Q <-Q
       model$R <-R
      }
    if(!is.null(model$z0)) model$z0 <- z
    if(!is.null(model$P0)) model$P0 <- P
    if(!is.null(model$rootP0)) model$rootP0 <- rP
    model
} #end setArrays.SS

setArrays.ARMA <- function(model, coefficients=NULL, constants=NULL) { 
	# N.B. Dimension and class info. is assumed accurate
       if (is.null(coefficients)) coefficients    <- coef(model)
                        else   model$coefficients <- coefficients
       a.pos  <- model$location
       i.pos  <- model$i
       j.pos  <- model$j
       const  <-if (is.null(constants)) model$const else constants #untested
       ca.pos <- model$const.location
       ci.pos <- model$const.i
       cj.pos <- model$const.j  
       l.pos  <- model$l
       cl.pos <- model$const.l
       A  <-  array(0,dim(model$A))
       B  <-  array(0,dim(model$B))
       m <- dim(model$C)[3]
       p <- dim(model$A)[2]
       TREND <- if(is.null(model$TREND)) NULL else rep(0,length(model$TREND))
       if(!is.null(m)) C  <-  array(0,dim(model$C))
       if (length(coefficients)>0) 
          {i <- a.pos == "A"
           A[cbind(l.pos[i],i.pos[i],j.pos[i])] <- coefficients[i]
           i <- a.pos == "B"
           B[cbind(l.pos[i],i.pos[i],j.pos[i])] <- coefficients[i]
           if(!is.null(m))
             {i <- a.pos == "C"
              C[cbind(l.pos[i],i.pos[i],j.pos[i])] <- coefficients[i]
             }
           i <- a.pos == "t"
           if(!is.null(TREND)) TREND[i.pos[i]] <- coefficients[i]
          }
    if (length(const)>0) 
          {i <- ca.pos == "A"
           A[cbind(cl.pos[i],ci.pos[i],cj.pos[i])] <- const[i]
           i <- ca.pos == "B"
           B[cbind(cl.pos[i],ci.pos[i],cj.pos[i])] <- const[i]
           if(!is.null(m))
             {i <- ca.pos == "C"
              C[cbind(cl.pos[i],ci.pos[i],cj.pos[i])] <- const[i]
             }
           i <- ca.pos == "t"
           if(!is.null(TREND)) TREND[ci.pos[i]] <- const[i]
          }
      model$A <- A
      model$B <- B
      if(!is.null(m)) model$C <- C 
      model$TREND <- TREND
      model 
} #end setArrays.ARMA


############################################################

# following is a workaround for ar in ts library

DSE.ar <- function(data, ...) {
  #fix for ar in R ts library (so that univariate case also gives array result)
  # before R 1.9.0 required ts not stats
  res <- stats::ar(ts(data), ...)
  if (! is.array(res$ar)) res$ar <- array(res$ar, c(length(res$ar),1,1))
  res
  }

printTestValue <- function (x, digits = 16){
    cat("c( ")
    if (all(is.na(x)))       cat("NAs")
    else if (is.null(x))     cat("NULL")
    else if (is.logical(x))  cat(x)
    else if (!is.R())        print(x, digits = digits)
    else if (is.matrix(x)) {
       for (i in 1:nrow(x)) {
	 cat("\n      ")
         for (j in 1:ncol(x)) 
	   cat(formatC(x[i,j], digits=digits, format="g"), ", ")
	 }
       cat("), ", ncol(x), ", ", nrow(x), ")\n")
       }
    else for (i in 1:length(x)) cat(", ", formatC(x[i], digits=digits, format="g"))
    cat(")\n")
    invisible()
    }


############################################################

#  Model simulation functions (to generate data)   <<<<<<<<<<

############################################################



simulate <- function(model, ...)UseMethod("simulate")

simulate.TSestModel <- function(model, input=inputData(model),
			sd=NULL, Cov=NULL, ...)
  {if (is.null(sd) & is.null(Cov)) Cov <- model$estimates$cov
   simulate(TSmodel(model), input=input, sd=sd, Cov=Cov, ...)
  }

simulate.SS <- function(model, input=NULL,
                 start=NULL, freq=NULL, sampleT=100, 
                 noise=NULL, sd=1, Cov=NULL, rng=NULL, 
                 compiled=.DSEflags()$COMPILED, ...)
{#  (... further arguments, currently disregarded)
if (is.TSestModel(model)) model <- TSmodel(model)
if (!is.SS(model)) stop("simulate.SS expecting an SS TSmodel.")
if (!checkConsistentDimensions(model)) stop("The SS model is not correct.")

 FF<-    model$F
 G <-    model$G
 H <-    model$H
 m <- dim(G)[2]
 if(is.null(m)) m <-0
 n <- dim(FF)[2]
 p <- dim(H)[1]
 if (m!=0)
   {if( is.null(input)) stop("input series must be supplied for this model.")
    if (sampleT != Tobs(input) ) input <- tfTruncate(input, end=sampleT)
   }
 if (is.innov.SS(model))
   {K <-    model$K}
 else
   {Q <-    model$Q
    if (ncol(Q)<n) Q <- cbind(Q,matrix(0,n,n-ncol(Q))) # Q assumed square
    R <-    model$R
   }
 
 if(is.null(rng)) rng <- setRNG() # returns setting so don't skip if NULL
 else        {old.rng <- setRNG(rng);  on.exit(setRNG(old.rng))  }
 
if (!is.null(noise) && is.matrix(noise)) noise <- list(w=noise)

set.ts <- TRUE             
if (!is.null(start))
  {if (!is.null(freq))   tf <- list(start=start, frequency=freq)
   else
      {warning("start is specified but not freq. Using freq=1.") 
       tf <- list(start=start, frequency=1)
      }
  }  
else if( (!is.null(input))   && is.tframed(input))   tf <- tframe(input)
else if ((!is.null(noise$w)) && is.tframed(noise$w)) tf <- tframe(noise$w)
else set.ts <-  FALSE



# It would be better to use makeTSnoise here (lag=1 and w0<-c(noise$w0) ) and
#  possibly two calls to get e for non-innov models. However, this will affect
#  historical comparisons so it needs to be done carefully!
 e <-NULL
 if (is.null(noise))
   {if (!is.innov.SS(model)) 
      {w0 <- rnorm(n)
       e <- matrix(rnorm(sampleT*n),sampleT,n)
       w <- matrix(rnorm(sampleT*p),sampleT,p)
      }
    else 
      {w <- matrix(NA,sampleT,p)
       w0 <- rep(NA,p)
       if (!is.null(Cov))
         {if(length(Cov) == 1) Cov <- diag(Cov, p)
	  W <- t(chol(Cov))
	  w.help <- t(W %*% matrix(rnorm((sampleT+1)*p),nrow=p,ncol=sampleT+1))
	  w0 <- w.help[1,]
	  w <- w.help[-1,]
         }
       else
         {if (length(sd)==1) sd <-rep(sd,p)
          for (i in 1:p)
            {w0[i] <- rnorm(1,sd=sd[i])
             w[,i] <- rnorm(sampleT,sd=sd[i])
            }
         }
      }
   }
 else
   {#  noise is not null
   # already done if (is.matrix(noise)) noise <- list(w=noise)
   if (is.null(noise$w))
       stop("supplied noise structure is not correct. w must be specified.")
   if (is.null(noise$w0)) noise$w0 <- rep(0,p)
   if ((!is.innov.SS(model)) &&  is.null(noise$e))
       stop("supplied noise structure is not correct. e must be specified for non-innovations form models.")

    w0 <- noise$w0
    if (is.matrix(w0)) w0 <- rep(0,p) # see note in man re VAR compatiblity
    w <- noise$w
    e <- noise$e
    sampleT <- Tobs(w)
   }

 y <- matrix(0,sampleT,p)
 state <- matrix(0,sampleT,n)
 if(is.null(model$z0)) z <- rep(0,n) # initial state
 else  z <- model$z0        

 if (is.innov.SS(model)) 
   {z <-  c(FF%*% z) + c(K %*% w0)
    if (m !=0) z <-  z + c(G%*%input[1,]) 
    y[1,]  <-  c(H %*% z) + c(w[1,])
   }   
 else      
   {z <-  c(FF%*% z) + c(Q %*% w0)
    if (m !=0) z <-  z + c(G %*% input[1,])  
    y[1,]  <-  c(H %*% z) + c(R %*% w[1,])
   }
 state[1,]<-z
 # Note: the first period is done above for both compiled and S versions so
 #    that initial conditions do not have to be passed to compiled code.
 if (compiled)
   {if (is.innov.SS(model))
      {Q <-matrix(0,n,n)
       R <-matrix(0,p,p)
       e <- matrix(0,sampleT,n)
      }
    else K <- matrix(0,n,p)
    if (m==0) 
      {input <- matrix(0,sampleT,1)
       G <- matrix(0,n,1)
      }
 #   storage.mode(y)     <- "double"
 #   storage.mode(state) <- "double"
 #   storage.mode(input) <- "double"
 #   storage.mode(w) <- "double"
 #   storage.mode(e) <- "double"
 #   storage.mode(FF) <- "double"
 #   storage.mode(G) <- "double"
 #   storage.mode(H) <- "double"
 #   storage.mode(K) <- "double"
 #   storage.mode(Q) <- "double"
 #   storage.mode(R) <- "double"
   #r <- .Fortran(simss_sym, 
   r <- .Fortran("simss", 
                         y=if(is.double(y)) y else as.double(y), 
                         state=if(is.double(state)) state else as.double(state), 
                         as.integer(m),
                         as.integer(n),
                         as.integer(p), 
                         as.integer(sampleT),  
                         if(is.double(input)) input else as.double(input),
                         if(is.double(w)) w else as.double(w),
                         if(is.double(e)) e else as.double(e),
                         if(is.double(FF)) FF else as.double(FF),
                         if(is.double(G)) G else as.double(G),   
                         if(is.double(H)) H else as.double(H),
                         if(is.double(K)) K else as.double(K), 
                         if(is.double(Q)) Q else as.double(Q),	 
                         if(is.double(R)) R else as.double(R),    
                         as.integer(is.innov.SS(model)), 
			 PACKAGE="dse"
			 )[c("y","state")]
    y <- r$y
    state <- r$state
    if (m==0) input <- NULL
   }
 else
   {for (Time in 2:sampleT)  
      {z <-  c(FF%*% z) 
       if (m !=0) z <-  z + c(G%*%input[Time,]) 
       if (is.innov.SS(model))  z <-  z + c(K%*%w[Time-1,])
       else       z <-  z + c(Q%*%e[Time-1,])
       state[Time,] <- z
       if (is.innov.SS(model)) y[Time,]  <-  c(H %*% z) + c(w[Time,])   
       else      y[Time,]  <-  c(H %*% z) + c(R %*% w[Time,])
   }  }
 if (set.ts)
   { y     <- tframed(y,     tf=tf, names=seriesNamesOutput(model)) 
     state <- tframed(state, tf=tf )
   }
 else  seriesNames(y) <- seriesNamesOutput(model)
 TSdata(list(input=input,output=y, state=state, version=version, 
   model=model, description="data generated by simulate.ss",
   noise=list(w0=w0,w=w, e=e, rng=rng, Cov=Cov, sd=sd)))
}

simulate.ARMA <- function(model, y0=NULL, input=NULL, input0=NULL,
                start=NULL, freq=NULL, sampleT=100,
                noise=NULL, sd=1, Cov=NULL,
                rng=NULL, noise.model=NULL, 
                compiled=.DSEflags()$COMPILED, ...)
{# see details in help("ARMA") and help("simulate.ARMA")

if (is.TSestModel(model)) model <- TSmodel(model)
if (!is.ARMA(model)) stop("simulate.ARMA expecting an ARMA TSmodel.")
if (!checkConsistentDimensions(model)) stop("The ARMA model is not correct.")

A<-    model$A
B <-    model$B
C <-    model$C
m <- dim(C)[3]
if (is.null(m)) m <-0
p <- dim(A)[2]
a <-dim(A)[1]
b <-dim(B)[1]
cc <- if (is.null(C)) 0 else  dim(C)[1]

TREND <- model$TREND
if (p == length(TREND)) TREND <- t(matrix(TREND, p, sampleT))

if (m!=0)
   {if( is.null(input)) stop("input series must be supplied for this model.")
    if (sampleT != Tobs(input) ) input <- tfTruncate(input, end=sampleT)
   }
 
if (p==1) invA0 <- matrix(1/A[1,,],1,1)
else      invA0 <- solve(A[1,,])
for (l in 1:a) A[l,,] <- invA0 %*% A[l,,]      # set A(0) = I      
for (l in 1:b) B[l,,] <- invA0 %*% B[l,,] 
if (m!=0) for (l in 1:dim(C)[1]) C[l,,] <- invA0 %*% C[l,,]  
if(!is.null(TREND)) TREND <- t(invA0 %*% t(TREND))

set.ts <- TRUE             

if (!is.null(noise)) {
   if (is.matrix(noise)) noise <- list(w=noise)
   if (is.null(noise$w0)) noise$w0 <-matrix(0,b,p)
   }

if (!is.null(start))
  {if (!is.null(freq))   tf <- list(start=start, frequency=freq)
   else
      {warning("start is specified but not freq. Using freq=1.") 
       tf <- list(start=start, frequency=1)
      }
  }  
else if( (!is.null(input))   && is.tframed(input))   tf <- tframe(input)
else if ((!is.null(noise$w)) && is.tframed(noise$w)) tf <- tframe(noise$w)
else set.ts <-  FALSE


noise <- makeTSnoise(sampleT,p,b, noise=noise, rng=rng,
                        Cov=Cov, sd=sd, 
			noise.model=noise.model,
                        start=start, frequency=freq)
if (is.null(sampleT)) sampleT<-noise$sampleT
 
 if(is.null(y0)) y0<-matrix(0,a,p)
 if((m!=0) & is.null(input0)) input0 <- matrix(0,dim(C)[1],m)
 y <- matrix(0,sampleT,p)  
 if (compiled)
   {if (m==0) 
      {input <- matrix(0,sampleT,1)
       input0 <- matrix(0,1,1)
       C <- matrix(0,1,1)
      }
    if (is.null(TREND)) TREND<- matrix(0,sampleT, p)
#    yo<- list(y=y, y0,m,p, a, b, cc, sampleT,input,input0,w,w0,A,B, C,TREND)
#    storage.mode(y)     <- "double"
#    storage.mode(y0)     <- "double"
#    storage.mode(input0)     <- "double"
#    storage.mode(noise$w)     <- "double"
#    storage.mode(noise$w0)     <- "double"
#    storage.mode(A)     <- "double"
#    storage.mode(B)     <- "double"
#    storage.mode(C)     <- "double"
#    storage.mode(TREND)     <- "double"

    #y <- .Fortran(simrma_sym, 
     y <- .Fortran("simrma", 
                        y=if(is.double(y)) y else as.double(y), 
                         if(is.double(y0)) y0 else as.double(y0),
                         as.integer(m),
                         as.integer(p), 
                         as.integer(a), 
                         as.integer(b), 
                         as.integer(cc), 
                         as.integer(sampleT),  
                         as.double(input[1:sampleT,]),
                         if(is.double(input0)) input0 else as.double(input0),
                         if(is.double(noise$w)) noise$w else as.double(noise$w),
                         if(is.double(noise$w0)) noise$w0 else as.double(noise$w0),
                         if(is.double(A)) A else as.double(A),
                         if(is.double(B)) B else as.double(B),   
                         if(is.double(C)) C else as.double(C),
                         if(is.double(TREND)) TREND else as.double(TREND), 
			 PACKAGE="dse"
			 ) [["y"]]

    if (m==0) 
      {input  <- NULL
       input0 <- NULL
      }
   }
 else
  {w0 <- noise$w0
   w<-noise$w
   if(!is.null(TREND)) y <- TREND  
   for (Time in 1:sampleT)  
   {for (l in 2:a) 
       if(Time+1-l<=0)
          if (p==1) y[Time,] <- y[Time,]-c(A[l,,]  *  y0[l-Time,]) 
          else      y[Time,] <- y[Time,]-c(A[l,,] %*% y0[l-Time,])
       else                
          if (p==1) y[Time,] <- y[Time,]-c(A[l,,]  *  y[Time+1-l,])
          else      y[Time,] <- y[Time,]-c(A[l,,] %*% y[Time+1-l,])
    for (l in 1:b) 
       if (Time+1-l<=0)
          if (p==1) y[Time,]<- y[Time,] +c(B[l,,]  *  w0[l-Time,]) 
          else      y[Time,]<- y[Time,] +c(B[l,,] %*% w0[l-Time,])
       else
          if (p==1) y[Time,]<- y[Time,] +c(B[l,,]  *  w[Time+1-l,])
          else      y[Time,]<- y[Time,] +c(B[l,,] %*% w[Time+1-l,])
    if (m!=0) for (l in 1:cc) 
       if (Time+1-l<=0)
          if (m==1) y[Time,]<- y[Time,] + c(C[l,,]  *  input0[l-Time,])
          else      y[Time,]<- y[Time,] + c(C[l,,] %*% input0[l-Time,])
       else
          if (m==1) y[Time,]<-y[Time,] + c(C[l,,]  *  input[Time+1-l,])
          else      y[Time,]<-y[Time,] + c(C[l,,] %*% input[Time+1-l,])
  }}
 if (set.ts) y <- tframed(y, tf=tf, names=seriesNamesOutput(model)) 
 else seriesNames(y) <- seriesNamesOutput(model)
 TSdata(list(input=input,output=y, 
          model=model, input0=input0, 
          description="data generated by simulate.ARMA", 
#          prior.args=prior.args, post.args=post.args,
          noise=noise))
}

#######################################################################

#functions which work on data (and models)  <<<<<<<<<<

############################################################

#     likelihood and residual calculation functions  <<<<<<<<<<

############################################################

residualStats <- function(pred, data, sampleT=nrow(pred), warn=TRUE){  
   e <- if (is.null(pred))     -data[1:sampleT,,drop=FALSE]
        else if (is.null(data)) pred[1:sampleT,,drop=FALSE]
        else               pred[1:sampleT,,drop=FALSE] - data[1:sampleT,,drop=FALSE]
   p <- ncol(e)
   #Om <-t(e) %*% e /sampleT
   Om <-crossprod(e)/sampleT
   if (any(is.na(Om))) {like1 <- like2 <- 1e100}  
   else if (any(Om >1e100)) {like1 <- like2 <- 1e100}  
   else
     {v <- La.svd(Om) #eigenvalues are not robust to degenerate density.
#     i <- v$d!=0
     i <- v$d > (v$d[1]*sqrt(.Machine$double.eps))
      if (!all(i))
        {if(warn) warning("The cov. matrix is singular. Working on subspace.")
         v$d <- v$d[i]
      #   v$u <- v$u[i,i, drop=FALSE]
      #   v$vt <- v$vt[i,i, drop=FALSE]
      #   e <- e[,i, drop=FALSE]
        }
      like1 <- 0.5 * sampleT * log(prod(v$d)) # det
      # $vt is transposed in La.svd, but not v in svd
#      if (1 == length(v$d)) OmInv <-  v$u %*% (1/v$d) %*% v$vt 
#      else OmInv <-  v$u %*% diag(1/v$d) %*% v$vt # more robust than solve(Om)
#  following is equivalent
#      OmInv <-  v$u %*% sweep(v$vt,1,1/v$d, "*") 
#      like2 <- sum(e * (e %*% OmInv)) /2
#  but this works out to (sampleT*p/2) and fixing for degenerate distributions:
       like2 <- (sampleT*length(v$d))/2
     }
   const <- (sampleT * p * log(2 * pi))/2
   invisible(list(like=c(const+like1+like2, const, like1,like2),
                  cov =Om, pred=pred, sampleT=sampleT))
   }

sumSqerror <- function(coefficients, model=NULL, data=NULL, error.weights=NULL) {
 if ( is.null(model)) stop("model missing") 
 if ( is.null(data))  stop("data missing") 
 if ( is.null(error.weights)) stop("error.weights missing")  
 sum(l(setArrays(model,coefficients=coefficients), data,
       result="weighted.sqerror",error.weights=error.weights))
 }

l <- function(obj1, obj2, ...)UseMethod("l")
l.TSdata <- function(obj1, obj2, ...) {l(obj2, obj1, ...) }
l.TSestModel <- function(obj1, obj2, ...) {l(TSmodel(obj1), obj2, ...)}


l.ARMA <- function(obj1, obj2, sampleT=NULL, predictT=NULL,result=NULL,
                error.weights=0,  compiled=.DSEflags()$COMPILED, 
		warn=TRUE, return.debug.info=FALSE, ...)
{#  (... further arguments, currently disregarded)
  #  see help("l.ARMA")
model <- if (is.TSestModel(obj1)) TSmodel(obj1) else obj1
if (!is.ARMA(model)) stop("l.ARMA expecting an ARMA TSmodel.")

#dat <- freeze(obj2)
dat <- obj2
tf <- tfTruncate(tframe(outputData(dat)), end=predictT)

if(!checkConsistentDimensions(model,dat)) stop("dimension error")
if (is.null(sampleT))  sampleT  <- TobsOutput(dat)
if (is.null(predictT)) predictT <- sampleT
if (sampleT > predictT) stop("sampleT cannot exceed predictT.")
if (sampleT > TobsOutput(dat)) stop("sampleT cannot exceed length of data.")

if (0 != nseriesInput(dat))
  {if (TobsInput(dat) < predictT)
      stop("input data must be at least as long as requested prediction.")
   if (any(is.na(inputData(dat)))) stop("input data cannot contain NAs.")
   if(!all(startInput(dat) == startOutput(dat)))
	   stop("input and output data must start at the same time.")
  }
if (any(is.na(outputData(dat)))) stop("output data cannot contain NAs.")

u <- inputData(dat)
y <- outputData(dat)
A<-    model$A
B <-   model$B
C <-   model$C
m <- dim(C)[3]
if (is.null(m)) m <-0
p <- dim(A)[2]
a <-dim(A)[1]
b <-dim(B)[1]

TREND <- model$TREND
if (p == length(TREND)) TREND <- t(matrix(TREND, p, predictT))

if (m == 0) 
   {if(!is.null(u)) 
      stop("model parameters indicate an no input but input data exists!")
   }

if (compiled)
  {if (m==0)
     {C <- array(0,c(1,p,1))    # can't pass 0 length array to compiled
      u <- matrix(0,predictT,1)
     }
   else # since input and output are declared same dim in fortran
      if (Tobs(y) < Tobs(u)) y <- 
	   rbind(y, matrix(0, Tobs(u) - Tobs(y), nseries(y)))
   if (is.null(TREND)) TREND <- matrix(0,predictT, p)
   is  <- max(m,p)

#   storage.mode(error.weights)     <- "double"
#   storage.mode(u)     <- "double"
#   storage.mode(y)     <- "double"
#   storage.mode(A)     <- "double"
#   storage.mode(B)     <- "double"
#   storage.mode(C)     <- "double"
#   storage.mode(TREND)     <- "double"

   #r  <- .Fortran(arma_sym,
   r  <- .Fortran("arma",
                         pred=matrix(1e20,predictT,p), # pred,     
                         as.integer(length(error.weights)),
                         weighted.sqerror=matrix(0,sampleT,p),
                         error.weights= if(is.double(error.weights)) 
			           error.weights else as.double(error.weights),
                         as.integer( m), 
                         as.integer( p) ,      
                         as.integer( dim(A)[1]),  # 1+order of A  
                         as.integer( dim(B)[1]),  # 1+order of B  
                         as.integer( dim(C)[1]),  # 1+order of C  
                         sampleT=as.integer(sampleT),
                         predictT=as.integer(predictT),
                         as.integer(Tobs(y)), 
                         if(is.double(u)) u else as.double(u), # as.double() works ok with compiled but
                          #messes up the dim(u) returned in the list
                         if(is.double(y)) y else as.double(y),	     
                         if(is.double(A)) A else as.double(A),  
                         if(is.double(B)) B else as.double(B),   
                         if(is.double(C)) C else as.double(C),
                         if(is.double(TREND)) TREND else as.double(TREND),
                         as.integer(is),  # scratch array dim
                         matrix(double(1),is,is),  # scratch array
                         matrix(double(1),is,is),  # scratch array
                         double(is),         # scratch array
                         integer(is*is),         # scratch array IPIV
                         PACKAGE="dse"
			 ) [c("pred", "weighted.sqerror")]
   if (all(0==error.weights)) r$weighted.sqerror <- NULL
  }
else   # start S version
  {prederror <- matrix(0,sampleT,p)  
   wt.err <- NULL
   invB0 <- solve(B[1,,])
   for (l in 1:a) A[l,,] <- invB0 %*% A[l,,]      # set B(0) = I      
   for (l in 1:b) B[l,,] <- invB0 %*% B[l,,]  
   if (m!=0) for (l in 1:dim(C)[1]) C[l,,] <- invB0 %*% C[l,,]  
   if(!is.null(TREND)) TREND <- t(invB0 %*% t(TREND))
   if (1 < length(error.weights)) wt.err <- matrix(0,predictT,p)    
   for (Time in 1:sampleT)  
      {if(!is.null(TREND)) vt <- -TREND[Time,]
       else vt    <-  rep(0,p) 
       for (l in 1:a)
          if (l<=Time)  #this is cumbersome but drop=FALSE leaves A,B,C as 3 dim arrays
              if (p==1) vt <- vt + c(A[l,,]  *  y[Time+1-l,])
              else      vt <- vt + c(A[l,,] %*% y[Time+1-l,])  
       if (b >= 2) for (l in 2:b) 
          if (l<=Time) 
             if (p==1) vt <- vt - c(B[l,,]  *  prederror[Time+1-l,]) 
             else      vt <- vt - c(B[l,,] %*% prederror[Time+1-l,])
       if (m!=0) for (l in 1:dim(C)[1]) 
          if (l<=Time) 
             if (m==1) vt <- vt - c(C[l,,]  *  u[Time+1-l,])
             else      vt <- vt - c(C[l,,] %*% u[Time+1-l,])  
       prederror[Time,] <- vt #this is not really the pred error unless B0=I   

       if (any(0!=error.weights))          
        {wt.err[Time,] <- error.weights[1]*vt^2  # weighted sq prediction error
         if (length(error.weights)>1)
           {for (h in 2:length(error.weights))
            if ( (Time+h-1) <= sampleT)
              {if(!is.null(TREND)) vt <- -TREND[Time,]
               else vt    <-  rep(0,p) 
               for (l in 1:a)
                  if (l < Time+h) 
                     if (p==1) vt <- vt + c(A[l,,]  *  y[Time+h-l,])
                     else      vt <- vt + c(A[l,,] %*% y[Time+h-l,])  
               if (b >= 2) for (l in 2:b) 
                  if (l < Time+h) 
                     if (p==1) vt <- vt - c(B[l,,]  *  prederror[Time+h-l,]) 
                     else      vt <- vt - c(B[l,,] %*% prederror[Time+h-l,])
               if (m!=0) for (l in 1:dim(C)[1]) 
                   if (l < Time+h) 
                      if (m==1) vt <- vt - c(C[l,,]  *  u[Time+h-l,])
                      else      vt <- vt - c(C[l,,] %*% u[Time+h-l,]) 
               wt.err[Time,] <- wt.err[Time,] + error.weights[h]*(solve(invB0 )%*%vt)^2
     }  }  }
  }

   pred <- matrix(0,predictT,p)
   pred[1:sampleT,] <- y[1:sampleT,,drop=FALSE] - prederror[1:sampleT,,drop=FALSE] %*% t(solve(invB0)) 
   # now multi-step predictions to predictT
   if (predictT > sampleT)
    {for (Time in (sampleT+1):predictT)  
      {invA0 <- solve(A[1,,])
       for (l in 1:a) A[l,,] <- invA0 %*% A[l,,]      # set A(0) = I      
       for (l in 1:b) B[l,,] <- invA0 %*% B[l,,]  
       if (m!=0) for (l in 1:dim(C)[1]) C[l,,] <- invA0 %*% C[l,,]  
       if(!is.null(TREND)) TREND <- t(invA0 %*% t(TREND))
       if(!is.null(TREND)) pred[Time,] <- pred[Time,]+TREND[Time,]
       if (a >= 2) for (l in 2:a) 
          if(Time+1-l<=sampleT)
             if (p==1) pred[Time,] <- pred[Time,]-c(A[l,,]  *  y[Time+1-l,]) 
             else      pred[Time,] <- pred[Time,]-c(A[l,,] %*% y[Time+1-l,])
          else                
             if (p==1) pred[Time,] <- pred[Time,]-c(A[l,,]  *  pred[Time+1-l,])
             else      pred[Time,] <- pred[Time,]-c(A[l,,] %*% pred[Time+1-l,])
       if (b >= 2) for (l in 2:b) 
          if (Time+1-l<=sampleT)
             if (p==1)  pred[Time,] <- pred[Time,] +c(B[l,,]  *  prederror[Time+1-l,]) 
             else       pred[Time,] <- pred[Time,] +c(B[l,,] %*% prederror[Time+1-l,])
       if (m!=0) for (l in 1:dim(C)[1]) 
          if (m==1) pred[Time,] <- pred[Time,] + c(C[l,,]  *  u[Time+1-l,]) 
          else      pred[Time,] <- pred[Time,] + c(C[l,,] %*% u[Time+1-l,])
      }
    }
   r<-list(pred=pred,   weighted.sqerror=wt.err)
  } # end of S version

tframe(r$pred) <- tf
if (! is.null(r$weighted.sqerror))  tframe(r$weighted.sqerror) <- tf
if((!is.null(result)) && (result == "pred")) return(r$pred)
r <- append(residualStats(r$pred, y, sampleT, warn=warn), 
        list(error.weights=error.weights, weighted.sqerror=r$weighted.sqerror))

if (return.debug.info) 
    r$debug.info <-list(m=m, p=p, a=a,b=b,c=c, A=A,B=B,C=C,TREND=TREND, 
          u=u, y=y,
          prederror=prederror, pred=pred, invB0=invB0, wt.err=wt.err, 
          error.weights=error.weights, sampleT=sampleT)

if ( is.null(result)) return(classed( # TSestModel constructor (l.ARMA)
              list(estimates=r, data=dat, model=model), "TSestModel"))
else 
   {if (result =="like") return(r$like[1]) # neg.log.like. from residualStats
    else { return(r[[result]]) }
   }
stop("should never get to here in l.ARMA.")
}                       


l.SS <- function(obj1, obj2, sampleT=NULL, predictT=NULL, error.weights=0,
                 return.state=FALSE, return.track=FALSE, result=NULL, 
		 compiled=.DSEflags()$COMPILED,
                 warn=TRUE, return.debug.info=FALSE, ...)
{#  (... further arguments, currently disregarded)
 # see  help("l.SS") and help("SS")
# could check that Q is symmetric  or positive definite but ...

model <- if (is.TSestModel(obj1)) TSmodel(obj1) else obj1
if (!is.SS(model)) stop("l.SS expecting an SS TSmodel.")

#data <- freeze(obj2)
data <- obj2
if(!checkConsistentDimensions(model, data)) stop("dimension error.\n")
if (is.null(sampleT))  sampleT  <- TobsOutput(data)
if (is.null(predictT)) predictT <- sampleT
if (sampleT > predictT) stop("sampleT cannot exceed predictT.\n")
if (sampleT > TobsOutput(data))
    stop("sampleT cannot exceed length of data.\n")
if (0 != nseriesInput(data))
  {if (TobsInput(data) < predictT)
      stop("input data must be at least as long as requested prediction.\n")
   if (any(is.na(inputData(data)))) stop("input data cannot contain NAs.\n")
  }
if (any(is.na(outputData(data)))) stop("output data cannot contain NAs.\n")

Innov <- is.innov.SS(model)
if (Innov & return.track) 
   warning("Tracking error is zero for an innovations model. track will not be calculated.")

FF<-    model$F
H <-    model$H
n <- dim(FF)[2]
p <- dim(H)[1]
y <- outputData(data)
m <- if(is.null(model$G)) 0 else dim(model$G)[2]

if (m != 0)
  {G <-model$G
   u <- inputData(data)
  } 
if (Innov)            # K or Q,R can be NUll in model, which messes up compiled
   {K <-    model$K
    Q <-    matrix(0,1,1)      #not used
    R <-    matrix(0,1,1)      #not used
    track <-array(0,c(1,1,1))  #not used
   }
else
   {Q <-    model$Q
    if (ncol(Q)<n) Q <- cbind(Q,matrix(0,n,n-ncol(Q))) # Q assumed square in compiled
    R <-    model$R
    K <-    matrix(0,n,p)      # this is used
    if(return.track) track <-array(0,c(predictT,n,n))
    else             track <-array(0,c(1,1,1))  #not used
   }
if (return.state | return.debug.info) state <- matrix(double(1),predictT,n)
else              state <- matrix(double(1),1,1)    #not used
if(is.null(model$z0)) z <-rep(0,n)   # initial state
else  z <-model$z0

# initial state tracking error  --  not used in innov. models
 P <-  if(!is.null(model$rootP0)) t(model$rootP0) %*% model$rootP0
       else if(!is.null(model$P0))  model$P0
       else   diag(1,n)

if (compiled)
  {if (m == 0) {
      G<-matrix(0,n,1)       # can't call compiled with 0 length arrays
      u <- matrix(0,predictT,1)
      } 
   else {# since input and output are declared same dim in fortran
      if (Tobs(y) < Tobs(u)) y <- 
	   rbind(y, matrix(0, Tobs(u) - Tobs(y), nseries(y)))
      }
#   storage.mode(error.weights)     <- "double"
#   storage.mode(state) <- "double"
#   storage.mode(track) <- "double"
#   storage.mode(u)     <- "double"
#   storage.mode(outputData(data))     <- "double"
#   storage.mode(FF)     <- "double"
#   storage.mode(G)     <- "double"
#   storage.mode(H)     <- "double"
#   storage.mode(K)     <- "double"
#   storage.mode(Q)     <- "double"
#   storage.mode(R)     <- "double"
#   storage.mode(z)     <- "double"
#   storage.mode(P)     <- "double"
IS <- max(n,p)
   #r <- .Fortran(kf_sym,
   r <- .Fortran("kf",
                  pred=matrix(double(1),predictT,p), #note, as.double messes dim of pred    
                  as.integer(length(error.weights)), 
                  weighted.sqerror=matrix(0,sampleT,p),
                  error.weights=if(is.double(error.weights))
		                 error.weights else as.double(error.weights),   
                  as.integer(return.state),
                  state=state,         
                  as.integer(return.track & !Innov),
                  track=if(is.double(track)) track else as.double(track),                  
                  as.integer(m), 
                  as.integer(n), 
                  as.integer(p), 
                  sampleT=as.integer(sampleT), 
                  predictT=as.integer(predictT), 
                  as.integer(Tobs(y)),  
                  if(is.double(u)) u else as.double(u), 
                  if(is.double(y)) y else as.double(y),  
                  if(is.double(FF)) FF else as.double(FF),   
                  if(is.double(G)) G else as.double(G),	
                  if(is.double(H)) H else as.double(H),  
                  if(is.double(K)) K else as.double(K), 
                  if(is.double(Q)) Q else as.double(Q),	   
                  if(is.double(R)) R else as.double(R),	 
                  as.integer(Innov),
                  if(is.double(z)) z else as.double(z),
                  if(is.double(P)) P else as.double(P),
	          as.integer(IS),           # scratch arrays for KF, IS
	          matrix(double(1),IS,IS),  #A
	          matrix(double(1),IS,IS),  # AA
	          matrix(double(1),IS,IS),  # PP
	          matrix(double(1),n,n),  # QQ
	          matrix(double(1),p,p),  # RR 
	          rep(double(1),IS),  # Z
	          rep(double(1),IS), # ZZ
	          rep(double(1),IS), # WW		   
                  integer(IS*IS),         # scratch array IPIV
		  PACKAGE="dse"
		  ) [c("pred","state","track","weighted.sqerror")]
   if (all(0==error.weights)) r$weighted.sqerror <- NULL
  }
else                  #  S version
  {vt    <-  rep(0,p)     # initial prediction error
   pred  <- matrix(0,predictT,p) 
   wt.err <- NULL
   if (1 < length(error.weights)) wt.err <- matrix(0,predictT,p)

## This is the beginning of the  most important part of the algorithm ####

   if ( ! Innov ) 
     {RR <- R %*% t(R)  
      QQ <- Q %*% t(Q)  
     }                            
                                   
   for (Time in 1:sampleT)  {
       if ( ! Innov) 
         {PH  <-  P %*% t(H)
          ft    <- ( H %*% PH )  + RR         
          ft    <-  (ft + t(ft))/2   # force ft to be symmetric 
          K   <-  t(solve(ft,t(FF %*% PH)))  
          P   <-  (FF %*% P %*% t(FF) ) - ( K %*% H %*% P %*% t(FF) ) + QQ  # P(t|t-1)
          P   <-  (P + t(P))/2  # force symmetry (eliminate rounding error problems)
          if (return.track) track[Time,,] <- P   # P(t|t-1)
          # note P(t|t) = P-P%*%t(H)%*%solve(H%*%t(H)+RR)%*%P
         } 
         
       z <- c(FF%*%z) + c(K%*%vt)  # E[z(t)| t-1 ]
       if (m !=0) z <- z + c(G%*%u[Time,])
       if (return.state | return.debug.info) state[Time,] <- z
       pred[Time,] <- Ey  <-  c(H %*% z)       # predicted output     
       vt<-  y[Time,] - Ey                     # prediction error 

## This is the end of the most important part of the algorithm ####

       if (any(0!=error.weights))          
        {wt.err[Time,] <- error.weights[1]*vt^2  # weighted sq prediction error
         if (length(error.weights)>1)
          {zh <- z
           for (h in 2:length(error.weights))
            if ( (Time+h-1) <= sampleT)
              {zh <-  c(FF%*%zh)
               if (h==2)  zh <- zh + c(K%*%vt) # vt is 0 for h>2
               if (m !=0) zh <- zh + c(G%*%u[Time+h-1,])
               wt.err[Time,] <- wt.err[Time,] + 
                          error.weights[h]*(y[Time+h-1,] -  c(H %*% zh))^2
   }   }  }   }


#   prederror <- y[1:sampleT,,drop=FALSE]-pred[1:sampleT,,drop=FALSE]

   # now multi-step prediction to predictT 
   # This requires u but not y (y is ignored if it is supplied)
   if (predictT > sampleT)
    {for (Time in (sampleT+1):predictT)  
       {z <-  c(FF%*% z) 
        if (m !=0) z <-  z + c(G%*%u[Time,]) 
        if (Time==sampleT+1) z <- z + c(K%*%vt)
        if (return.state) state[Time,] <- z
        pred[Time,]  <-  c(H %*% z)                  # predicted output 
       }
     }
   r<- list(pred=pred, state=state, track=track, weighted.sqerror=wt.err) 
  }   # end of S version

tf <- tfTruncate(tframe(outputData(data)), end=predictT)
tframe(r$pred) <- tf
if (! is.null(r$weighted.sqerror))  tframe(r$weighted.sqerror) <- tf

filter <-NULL
if (return.track) filter$track <- if (Innov) NULL else tframed(r$track, tf)
if (return.state) filter$state <- tframed(r$state,tf)
 
if((!is.null(result)) && (result == "pred")) return(r$pred)
r <- append(residualStats(r$pred, outputData(data), sampleT, warn=warn), 
        list(error.weights=error.weights, weighted.sqerror=r$weighted.sqerror))

if (return.debug.info) 
   r$debug.info <- list(m=m, p=p, c=c, G=G,FF=FF,K=K,P=P, H=H, u=u,y=y,
        pred=pred, wt.err=wt.err, error.weights=error.weights, sampleT=sampleT)

if ( is.null(result)) return( 
    classed(list(estimates=r, data=data, model=model, filter=filter),
            "TSestModel")) # constructor (l.SS)
else if (result =="like") return(r$like[1]) # neg.log.like. from residualStats
else return(r[[result]]) 
  
"should never get to here in l.SS"
} # end of l.SS



smoother <- function(model, data, compiled=.DSEflags()$COMPILED) UseMethod("smoother")

smoother.TSestModel <- function(model, data=TSdata(model),
                                compiled=.DSEflags()$COMPILED){
    r <- smoother(TSmodel(model), data=data, compiled=compiled)
    r$estimates <- model$estimates # otherwise convergence info is lost
    r
    }
    
smoother.TSmodel <- function(model, data, compiled=.DSEflags()$COMPILED){
    if  (!is.nonInnov.SS(model))
       stop("smoothing needs an nonInnov state space model (class nonInnov.SS.TSmodel).")
    # should neot get to here. (smoother.nonInnov should have been called
    # instead, but ...)
    smoother.nonInnov(model, data, compiled=compiled)
    }

smoother.nonInnov <- function(model, data, compiled=.DSEflags()$COMPILED){
 # See help("smoother") and help("SS") for details of the model: 
 data <- TSdata(data)
 z <-l(model, data, return.state=TRUE,return.track=TRUE, compiled=compiled)
 filter    <- z$filter
 estimates <- z$estimates
 if  (is.null(filter$state) |  is.null(filter$track)) 
     filter <- l(model,data, return.state=TRUE,return.track=TRUE)$filter
 model     <- TSmodel(model)
 if  (!is.nonInnov.SS(model)) stop("smoother expecting an nonInnov.SS TSmodel.")
 if (is.null(model$G))
  {m<-0
   G<-matrix(0,dim(model$F)[2],1)   # can't call compiled with 0 length arrays
   u <- matrix(double(1),nrow(filter$state),1)
  }
 else
  {m <- dim(model$G)[2]
   G <-model$G
   u <- inputData(data)
   if (! is.double(u)) u <- array(as.double(u), dim(u))                  
  } 
 sampleT  <-min(nrow(u), nrow(filter$state), dim(filter$track)[1])
 QQ <- model$Q %*% t(model$Q)           
 RR <- model$R %*% t(model$R) 
 n <- dim(model$F)[2]          
 p <- dim(model$H)[1]
 if (compiled)
   {#storage.mode(filter$state)     <- "double"
    #storage.mode(filter$track)     <- "double"
    #storage.mode(u)     <- "double"
    #storage.mode(outputData(data))     <- "double"
    #storage.mode(model$F)     <- "double"
    #storage.mode(G)     <- "double"
    #storage.mode(model$H)     <- "double"
    #storage.mode(RR)     <- "double"
    
    #note, as.double messes dim needed in result, but need attributes removed
    IS <- max(n,p)

    #r <- .Fortran(smooth_sym, 
    r <- .Fortran("smooth", 
                         state=array(as.double(filter$state), dim(filter$state)),      
                         track=array(as.double(filter$track), dim(filter$track)),     
                         u,	       # input
                         if(is.double(outputData(data))) outputData(data) else 
			        as.double(outputData(data)),
			 as.integer(n),
                         as.integer(m),
                         as.integer(p),  
                         sampleT =as.integer(sampleT), 
                         as.double(model$F),   
                         as.double(G),	
                         as.double(model$H),  
                         as.double(RR),	 
                         as.integer(IS),           # scratch array dim IS
			 matrix(double(1),IS,IS),   # scratch array A
                         matrix(double(1),IS,IS),   # scratch array D
                         matrix(double(1),IS,IS),   # scratch array L
                         matrix(double(1),IS,IS),   # scratch array PT1
                         double(IS),         # scratch array ZT
                         integer(IS*IS),         # scratch array IPIV
                         #chk1=array(double(1), c(sampleT, IS,IS)),     
                         #chk2=array(double(1), c(sampleT, IS,IS)),     
                         PACKAGE="dse"
			 )[c("state","track")]
   }
 else   # S version
   {FF<-  model$F
    H <-  model$H
    #   filter$state is the one step ahead state estimate E{z(t)| y(t-1), u(t)}
    #   zt below is the filter state estimate E{z(t)| y(t), u(t+1)}
    #   filter$track is tracking error P(t|t-1)
    # A5-A7 just use the filter information , which was calculated 1:sampleT.
    # A8-A10 calculate smooth state and track, which is done in reverse sampleT:1 
    sm <- array(NA,dim(filter$state))  # smoother state estimate
    sm[sampleT,] <- filter$state[sampleT,]
    tr <- array(NA,dim(filter$track))  # smoother tracking error
    tr[sampleT,,] <- filter$track[sampleT,,]
    for (Time in (sampleT-1):1) 
      {#K <- filter$track[Time,,] %*% t(H) %*% 
       #       solve(H %*% filter$track[Time,,] %*% t(H) + RR) #(A5)
       K <- filter$track[Time,,] %*%  
              t(solve(t(H %*% filter$track[Time,,] %*% t(H) + RR), H)) #(A5)
       zt <- filter$state[Time,] + K %*% 
       		(outputData(data)[Time,] - H %*% filter$state[Time,]) #(A6)
       #if (m!=0) zt <- zt + G %*% u[Time,] 
       P <- filter$track[Time,,] - K %*% H %*% filter$track[Time,,] #P(t|t) (A7)
       P <- (P+t(P))/2 #force symmetry to avoid rounding problems
       #J <- P %*% t(FF) %*% solve(filter$track[Time+1,,])             #(A8) 
       J <- P %*%  t(solve(t(filter$track[Time+1,,]), FF))             #(A8) 
       if (m==0)                                          
          sm[Time,] <- zt + J %*% (sm[Time+1,] - FF %*% zt)              #(A9) 
       else 
          sm[Time,] <- zt + J %*% (sm[Time+1,] - FF %*% zt - G %*% u[Time+1,]) #check   
       P <- P + J %*% (tr[Time+1,,] - filter$track[Time+1,,]) %*% t(J)    #(A10)
       tr[Time,,]<- (P+t(P))/2   #force symmetry to avoid rounding problems
      }
     r <- list(state=sm, track=tr) 
    }  # end S version
 
  tf <- tframe(outputData(data))
  names <- seriesNames(filter$state) 
  classed(list(estimates=estimates, data=data, model=model, 
               filter=filter, 
	       smooth=list(state=tframed(r$state, tf, names=names),
	                   track=tframed(r$track, tf))),
       "TSestModel") # constructor (smoother)
} # end of smoother
  

state <- function(obj, smoother=FALSE, filter=!smoother) {
     if (!inherits(obj,"TSestModel"))
        stop("The argument needs to be a TSestModel from an SS model.")
     if(filter & smoother) stop("only one of filter and smoother can be specified.")
     if(filter) if(!is.null(obj$filter$state)) return(obj$filter$state)
     else return(l(TSmodel(obj), TSdata(obj), return.state=TRUE)$filter$state)
     if (!is.null(obj$smooth)) return(obj$smooth$state)
     if (!inherits(TSmodel(obj),"nonInnov"))
          stop("A smoother state cannot be calculated . The argument needs to be a TSestModel from an nonInnov SS model.")
     state(smoother(obj), smoother=TRUE)
     }
 
#trackingError <- function(obj, ... should do this as for state 



############################################################

#     parameter estimation functions   <<<<<<<<<<

############################################################


estVARXls <- function(data, subtract.means=FALSE, re.add.means=TRUE,
   standardize=FALSE, unstandardize=TRUE, max.lag=NULL, 
   trend=FALSE, lag.weight=1.0, warn=TRUE) 
{# Estimate VAR model with exogenous variable using lsfit(). 
 # Returns a TSestModel.
 # This is very similar to estimating Markov parameters (see estSSMittnik).
 # Residuals,etc, are calculated by evaluating the estimated model with ARMA.
 # Data should be of class TSdata.
 # lag.weight is an exponential weight applied to lags. It should be in (0,1].
   if (is.null(max.lag)) max.lag <- 6
   #data <- freeze(data)
   names <- seriesNames(data)
   m <-  nseriesInput(data)
   p <- nseriesOutput(data)
   if (0 == p) stop("estVARXls requires output data to estimate a model.")
   missing.data <- any(is.na(outputData(data)))
   if (0 != m) { # ensure same time window
	inputData(data) <- tfwindow(inputData(data),
              start=startOutput(data), end=endOutput(data), warn=FALSE)
	outputData(data) <- tfwindow(outputData(data),
              start=startInput(data), end=endInput(data), warn=FALSE)
	missing.data <- any(missing.data, is.na(inputData(data)))
   }
   N <- TobsOutput(data)
   if (standardize)
     {svd.cov <- La.svd(var(outputData(data)))
      scalefac <- svd.cov$u %*% diag(1/svd.cov$d^.5, ncol=p)
      data <- scale(data, scale=list(output=scalefac))
     }
   if (subtract.means)
    {if(m!=0)
       {input.means<-colMeans(inputData(data))
        inputData(data)<-inputData(data)-t(matrix( input.means, m,N))
       }
     output.means <- colMeans(outputData(data))
     outputData(data)  <- outputData(data) - t(matrix(output.means, p,N))
    }
      # shift input to give feedthrough in one period
   if (m != 0) {z <- cbind(inputData(data)[2:N,],outputData(data)[1:(N-1),])}
   else z <- outputData(data)

 # The matrix Past is blocks of data:
 #  [in | out-1 | in-1 | out-2 | ... | in-max.lag | out-max.lag-1 ]
 # so the coef. matrix M has a corresponding structure.
 
   Past <- matrix(NA,N-max.lag,(p+m)*(max.lag))
   for (i in 0:(max.lag-1)) 
      Past[,(1+(m+p)*i):((m+p)*(1+i))] <-z[(max.lag-i):(N-1-i),] /(lag.weight^i)
   fit <- lsfit(Past,outputData(data)[(max.lag+1):N,,drop=FALSE],intercept=trend)
   if(missing.data)
     fit.res <-TSdata(output=fit$residual) #used only in the case of missing data for scaling
   M <- t(fit$coef)
   # correct exogenous blocks for normalization:
   if (standardize && (m!=0)) 
     {Tinv <- diag(svd.cov$d^.5, ncol=p)%*%svd.cov$u
      for (i in 0:(max.lag-1)) 
         M[,(1+(m+p)*i):(m+(m+p)*i)] <- Tinv %*% M[,(1+(m+p)*i):(m+(m+p)*i)]
     }
   TREND <- NULL
   if (trend)
     {TREND <- M[,1,drop=FALSE]
      M<-M[,2:(dim(M)[2]),drop=FALSE]
     }
   if (subtract.means & re.add.means)
    {if(m!=0) 
       {inputData(data)<-inputData(data) + t(matrix( input.means, m,N))
        Past<-Past +
           t(matrix(c(input.means,output.means),(p+m)*(max.lag),N-max.lag))
       }
     else
        Past <-Past+t(matrix(output.means, (p+m)*(max.lag),N-max.lag))
     outputData(data)  <- outputData(data) + t(matrix(output.means, p,N))
     # and correct estimation for non-zero mean:
     M<-M+estVARXmean.correction(Past, 
                       outputData(data)[(max.lag+1):N,,drop=FALSE],M, warn=warn)
    }
   A <- array(NA, c(1 + max.lag, p, p))
   A[1,,] <-diag(1, p)
   for (i in 0:(max.lag-1)) 
          A[2+i,,] <-  -M[,(m+1+(m+p)*i):(m+p+(m+p)*i),drop=FALSE] %*% 
                                 diag(lag.weight^i,p)
   if (m==0) C <- NULL   # no exog. variable
   else               # NB. there is an implicit shift in inputData(data)
      {C <-array(NA,c(max.lag,p,m))
       for (i in 0:(max.lag-1)) 
         C[1+i,,] <- M[,(1+(m+p)*i):(m+(m+p)*i),drop=FALSE]  %*% 
                                 diag(lag.weight^i,m)
      }
   B <- array(diag(1,p),c(1,p,p))
   model <-ARMA( description="model estimated by estVARXls",
              A=A,B=B,C=C,TREND=TREND)
   seriesNames(model) <- seriesNames(data)

   if (standardize & unstandardize)
     {scalefac <- solve(scalefac)
      data <- scale(data,  scale=list(output=scalefac))
      model<- scale(model, scale=list(output=scalefac))
      if(missing.data) fit.res <- scale(fit.res, scale=list(output=scalefac))
     }

   if(missing.data)
       {if (warn) warning(
          "missing data kludge. Predictions are reconstructed from lsfit residuals")
        return(fake.TSestModel.missing.data(model,data, fit.res$output,max.lag))
       }
   else return(l(model, data, warn=warn))
}

fake.TSestModel.missing.data <- function(model,data, residual, max.lag, warn=TRUE)
  {pred <- rbind(matrix(NA,max.lag,nseriesOutput(data)),residual) + 
                                                           outputData(data)
   r <- list(estimates = residualStats(pred, outputData(data), warn=warn),
               data = data, model = model)
   classed(r, "TSestModel") # fake missing data constructor
  }

estVARXmean.correction <- function(X, y, bbar,
                     fuzz=sqrt(.Machine$double.eps), warn=TRUE)
{# correction for model estimated with means subtracted
 Xbar <- t(array(colMeans(X), rev(dim(X))))
 ybar <- t(array(colMeans(y), rev(dim(y))))
 v <- La.svd(t(X)%*%X) # this is more robust than solve()
 if (warn && any(abs(v$d[1]*fuzz) > abs(v$d) ) ) 
   warning("The covariance matrix is nearly singular. Check for linearly related data.")
# if(1 == length(v$d))OmInv <- v$u %*% (1/v$d) %*% v$vt
# else OmInv <- v$u %*% diag(1/v$d) %*% v$vt	
#  following is equivalent
# OmInv <-  t(v$vt) %*% sweep(t(v$u),1,1/v$d, "*") CHECK
 OmInv <-  crossprod(v$vt, sweep(t(v$u),1,1/v$d, "*")) 

# t(-OmInv %*% ( 
#        (2*t(Xbar)%*%X - t(Xbar)%*%Xbar) %*% t(bbar) 
#         - t(Xbar)%*%y - t(X)%*%ybar + t(Xbar)%*%ybar ))
 t(-OmInv %*% ( 
        (2*crossprod(Xbar, X) - crossprod(Xbar)) %*% t(bbar) 
         - crossprod(Xbar, y) - crossprod(X,ybar) + crossprod(Xbar, ybar)))
}


estVARXar <- function(data, subtract.means=FALSE,  re.add.means=TRUE,
      standardize=FALSE, unstandardize=TRUE, aic=TRUE, max.lag=NULL, 
      method="yule-walker", warn=TRUE) 
{
   #data <- freeze(data)
   m <-  nseriesInput(data)
   p <- nseriesOutput(data)
   if (0 == p) stop("estVARXar requires output data to estimate a model.")
   N <- TobsOutput(data)
   if (standardize)
     {svd.cov <- La.svd(var(outputData(data)))
      scalefac <- svd.cov$u %*% diag(1/svd.cov$d^.5, ncol=p)
      data <- scale(data, scale=list(output=scalefac))
     }
   if(m!=0) input.means <- colMeans(inputData(data))
   output.means <- colMeans(outputData(data))
   if (subtract.means)
    {if(m!=0) inputData(data) <- inputData(data)-t(matrix( input.means, m,N))
     outputData(data)  <- outputData(data) - t(matrix(output.means, p,N))
    }
   if (m==0)  zdata <- outputData(data)  # no exog. variable
   else       zdata <- cbind(inputData(data),outputData(data))
         # NB. there is an implicit shift in inputData(data) in the line above
         #   because ar estimates lag parameters,
         # so C[1,,] is for one lag in u.
   if (is.null(max.lag))  AC<- DSE.ar(zdata, aic=aic, method=method)
   else                   AC<- DSE.ar(zdata, aic=aic, method=method, order.max=max.lag)
   if (re.add.means)
     {if (subtract.means)
        {if(m!=0) inputData(data)<-inputData(data) + t(matrix( input.means, m,N))
         outputData(data)  <- outputData(data) + t(matrix(output.means, p,N))
        }
      # and correct estimation for non-zero mean:
      max.lag <- AC$order
      if (max.lag==0) stop("all lags eliminated by AIC order selection.")
      if (m==0)  zdata <- outputData(data)  # no exog. variable
      else       zdata <- cbind(inputData(data),outputData(data))
      Past <- matrix(NA,N-max.lag,(p+m)*(max.lag))
      for (i in 0:(max.lag-1)) 
         Past[,(1+(m+p)*i):((m+p)*(1+i))] <-zdata[(max.lag-i):(N-1-i),]
      M <-estVARXmean.correction(Past, zdata[(max.lag+1):N,,drop=FALSE],
                  matrix(aperm(AC$ar, c(2,3,1)),m+p,(m+p)*max.lag), warn=warn)
      AC$ar<-AC$ar + aperm(array(M,c(m+p,m+p,max.lag)), c(3,1,2))
     }
   A <- array(0, c(1 + AC$order, p, p))
   A[1,,] <-diag(1, p)
   if (0==AC$order)
      warning("lagged output variables eliminated by AIC order selection.")
   else 
      A[2:(1+AC$order),,] <-  -AC$ar[, (m+1):(m+p), (m+1):(m+p), drop=FALSE]
   if (m==0) C <- NULL
   else 
     {if (0==AC$order)
             {warning("input variables eliminated by AIC order selection.")
              C <-array(0, c(1,p,m))
             }
      else C <-array(AC$ar[,(m+1):(m+p),1:m],c(AC$order,p,m))
     }
   B <- array(diag(1,p),c(1,p,p))
   model <-ARMA(
            description="model estimated by estVARXar",
            A=A,B=B,C=C, 
            names = seriesNames(data) )
 
   if (standardize & unstandardize)
     {scalefac <- solve(scalefac)
      data <- scale(data,  scale=list(output=scalefac))
      model<- scale(model, scale=list(output=scalefac))
     }
   l(model, data, warn=warn)
}

estSSfromVARX <- function(data, warn=TRUE, ...){
	# estimate a VARX model and convert to state space
	#  estimate a nested-balanced state space model by svd a la Mittnik from
	#   least squares estimate of VAR coefficents.
	model <-estVARXls(data, warn=warn, ...)
	l(toSS(model$model), model$data, warn=warn)
}



############################################################################

#  functions for model estimation (see also VARX ) and reduction   <<<<<<<<<<<<<

############################################################################


estWtVariables <- function(data, variable.weights,
                        estimation="estVARXls", estimation.args=NULL)
{ if (is.matrix(variable.weights))
    {if (any(La.svd(variable.weights)$d == 0))  
       stop("variable.weights transformation must be non singular.")
    }
 else
   {if (any(variable.weights == 0))  stop("variable.weights must be non zero.")
    variable.weights <- diag(variable.weights)
   }
 inv.wts <- solve(variable.weights)
 dimnames(inv.wts)          <-list(NULL, seriesNamesOutput(data))
 dimnames(variable.weights) <-list(NULL, seriesNamesOutput(data))
 scaled.model <- do.call(estimation, append(list(
                  scale(data, scale=list(output=inv.wts))),  estimation.args))
#           freeze(scale(data, scale=list(output=inv.wts)))), estimation.args))
 model <-scale(TSmodel(scaled.model), scale=list(output=variable.weights))
 model$description <- 
    paste("model estimated by estWtVariables with", estimation)
 l(model, data)
}

estMaxLik <- function(obj1, obj2=NULL, ...) UseMethod("estMaxLik")

estMaxLik.TSestModel <- function(obj1, obj2=TSdata(obj1), ...) {
	 # if obj1 is result from a previous estMaxLik then the gradient
	 # hessian and other information should be extracted, but
	 estMaxLik(TSmodel(obj1), obj2, ...) }

estMaxLik.TSdata <- function(obj1, obj2, ...) 
	 estMaxLik(TSmodel(obj2), obj1, ...) 

estMaxLik.TSmodel <- function(obj1, obj2, 
	algorithm="optim",
	algorithm.args=list(method="BFGS", upper=Inf, lower=-Inf, hessian=TRUE),
	...)
{# maximum likelihood estimation
 # "nml" algorithm.args=list(hessian=T, iterlim=20, 
 #     dfunc=gradNumerical, line.search="nlm",ftol=1e-5, gtol=1e-3,)
 Shape <- obj1
 #data <- freeze(obj2)
 data <- obj2
 func.like <- function(coefficients, Shape,data)
      {l(setArrays(Shape, coefficients=coefficients), data, result="like",
         warn=FALSE) }

 if (algorithm=="optim")
    {results <- optim(coef(Shape), func.like, method=algorithm.args$method,
	gr=algorithm.args$gr, 
	lower=algorithm.args$lower, upper=algorithm.args$upper,
	control=algorithm.args$control, hessian=algorithm.args$hessian,
	Shape, data) 
     parms <- results$par
     converged <- results$convergence == 0
     convergenceCode <- results$convergence
     description <- paste("Estimated with max.like/optim")
    }
 else if (algorithm=="nlm")
   {warning("This has not been tested recently (and there have been changes which may affect it.")
    results <-nlm(func.like, coef(Shape), hessian=algorithm.args$hessian, 
    	iterlim=algorithm.args$iterlim)
    parms <- results$estimate
    converged <- results$code <= 2
    convergenceCode <- results$code
    description <- paste("Estimated with max.like/nlm")
   }
  else stop(paste("Minimization method ", algorithm, " not supported."))
# else if (algorithm=="nlmin")
#   {warning("This has not been tested recently (and there have been changes which may affect it.")
#     results <-nlmin(func.like, coef(Shape), max.iter=algorithm.args$max.iter, 
#     	max.fcal=5*algorithm.args$max.iter, ckfc=0.01)
#     parms <- results$parms
#     converged <- results$converged
#     convergenceCode <- results$converged
#     # emodel$est$converged should be improved with conv.type info
#     description <- paste("Estimated with max.like/nlmin")
#   }
# else if (algorithm=="dfpMin")
#    {stop("This optimization method is no longer supported.")
#     results <- dfpMin(func.like, coef(Shape), 
#	dfunc=algorithm.args$dfunc, 
#	max.iter=algorithm.args$max.iter, 
#	ftol=algorithm.args$ftol, gtol=algorithm.args$gtol, 
#	line.search=algorithm.args$line.search) 
#     parms <- results$parms
#     converged <- results$converged
#     convergenceCode <- results$converged
#     description <- paste("Estimated with max.like/dfpMin")
#    }
 emodel <- l(setTSmodelParameters(setArrays(Shape, coefficients=parms)),data)
 emodel$estimates$algorithm <- algorithm
 emodel$estimates$results   <- results
 emodel$estimates$converged <- converged
 emodel$estimates$convergenceCode <- convergenceCode
 emodel$model$description <- paste(description, 
       if(!converged) " (not " else " (", "converged",
       ") from initial model: ", Shape$description)
 emodel
}    


estBlackBox <- function(data,...)
{# call current best black box technique.
  estBlackBox4(data, ...)
}


estBlackBox1 <- function(data,estimation="estVARXls", 
     reduction="MittnikReduction", criterion="taic", trend=FALSE, 
     subtract.means=FALSE, verbose=TRUE, max.lag=6)
{if ((estimation!="estVARXls") && (trend) )
     {cat("Trend estimation only support with estVARXls.\n")
      cat("Proceeding using estVARXls.\n")
      estimation<-"estVARXls"
     }

 if(estimation=="estVARXls")
     model <- estVARXls(data,trend=trend, subtract.means=subtract.means, 
                          max.lag=max.lag)
 else if(estimation=="estVARXar")
     model <- estVARXar(data, subtract.means=subtract.means, max.lag=max.lag)
 else if(estimation=="estVARXls")
     model <- estVARXls(data,trend=trend, subtract.means=subtract.means, 
                        max.lag=max.lag)
 else if(estimation=="estSSMittnik")
     model <- estSSMittnik(data,max.lag=max.lag, 
                             subtract.means=subtract.means, normalize=FALSE)
 else
   stop("estimation technique not supported.")
 if (verbose) 
   cat("First VAR model,              lags= ", dim(model$model$A)[1]-1,
       ", -log likelihood = ", model$estimates$like[1], "\n")
 model <- l(toSS(model),data)
 n <- dim(model$model$F)[1]
 if (verbose) cat("Equivalent    state space model, n= ", n,
                  ", -log likelihood = ", model$estimates$like[1], "\n")
 if (1 < n)
   {model <- do.call(reduction,
                     list(model, criterion=criterion, verbose=verbose))
   #model <- eval(call(reduction,model,criterion=criterion, verbose=verbose))
    if (verbose) 
       cat("Final reduced state space model, n= ", nstates(model),
           ", -log likelihood = ", model$estimates$like[1], "\n")
   }
  if (verbose &&  dev.cur() != 1 ) checkResiduals(model)
 model
}


estSSMittnik <- function(data, max.lag=6, n=NULL, 
    subtract.means=FALSE, normalize=FALSE)
{ #data <- freeze(data)
  m <- ncol(inputData(data))
  if(is.null(m))  m <- 0
  p <- ncol(outputData(data))
  if (0 == p) stop("estSSMittnik requires output data to estimate a model.")
  N <- nrow(outputData(data))
  if (subtract.means)
    {if(m!=0)inputData(data)<-inputData(data)-t(matrix(colMeans(inputData(data)), m,N))
     outputData(data)<- outputData(data) - t(matrix(colMeans(outputData(data)), p,N))
    }
  if (normalize)
    {svd.cov <- La.svd(var(outputData(data)))
     outputData(data) <- outputData(data) %*% svd.cov$u %*% diag(1/svd.cov$d^.5)
    }
      # shift input to give feedthrough in one period
  if (m != 0) {z <- cbind(inputData(data)[2:N,],outputData(data)[1:(N-1),])}
  else z <- outputData(data)
  Past <- matrix(NA,N-1-max.lag,(p+m)*(1+max.lag))
  for (i in 0:max.lag) 
    Past[,(1+(m+p)*i):((m+p)*(1+i))] <-z[(1+max.lag-i):(N-1-i),]
  M <- t(lsfit(Past,outputData(data)[(max.lag+2):N,],intercept=FALSE)$coef)
  if (normalize && (m!=0))  # correct exogenous blocks for normalization
    {Tinv <- diag(svd.cov$d^.5)%*%svd.cov$u
     for (i in 0:max.lag) 
       M[,(1+(m+p)*i):(m+(m+p)*i)] <- Tinv %*% M[,(1+(m+p)*i):(m+(m+p)*i),drop=FALSE]
    }
  if (p==1) M <-matrix(M,1,length(M))
#  browser()
  model <-SVDbalanceMittnik( M, m=m, n=n )$model
  z <-"nested-balanced model from least sq. estimates of markov parameters a la Mittnik"
  if(subtract.means) z <-paste(z," - means subtracted")
  if(normalize)      z <-paste(z," - outputs normalized")
  model$description  <-z
  seriesNames(model)  <- seriesNames(data)
  l(model, data)
}


MittnikReduction <- function(model, data=NULL, criterion=NULL, 
      verbose=TRUE,warn=TRUE)
{if (is.TSestModel(model)) 
    {if (is.null(data)) data <-TSdata(model)
     model <- TSmodel(model)
    }
  if(!is.TSmodel(model)) stop("MittnikReduction expecting a TSmodel.")
  if(is.null(data))
    stop("Reduction requires data to calculate criteria (balancing does not).")
  nMax <- if(is.SS(model)) dim(model$F)[2] else NULL
  MittnikReduction.from.Hankel( markovParms(model), nMax=nMax, data=data, 
        criterion=criterion, verbose=verbose, warn=warn)
}

MittnikReduction.from.Hankel <- function(M, data=NULL, nMax=NULL, 
   criterion=NULL, verbose=TRUE, warn=TRUE) 
#   Spawn=if (exists(".SPAWN")) .SPAWN else FALSE)
{ # See the documentation for MittnikReduction.

   read.int <- function(prmt)
     {err <-TRUE
      while (err)
       {cat(prmt)
 	n <- as.integer(readline()) # crude. this truncates reals
 	if (is.na(n)) cat("value must be an integer\n")
 	else err <- FALSE
       }
      n
     }

   #data <- freeze(data)
   m <-ncol(inputData(data))      # dim of input series
   if(is.null(m))m<-0
   z <-SVDbalanceMittnik(M, m, nMax)
   largeModel <- z$model
   svd.crit    <-z$crit
   n <- dim(largeModel$F)[1]
#   if (!Spawn)
#     {
      values <- NULL 
      for (i in 1:n) 
        {if(m!=0) z <-largeModel$G[1:i,,drop=FALSE]
         else     z <-NULL
         z <-SS(F.=largeModel$F[1:i,1:i,drop=FALSE],G=z,
                  H=largeModel$H[,1:i,drop=FALSE],K= largeModel$K[1:i,,drop=FALSE])
         z <-informationTestsCalculations(l(z,data, warn=warn))
         values <-rbind(values, z)
         if (verbose) cat(".")
        }
 #     }
#   else 
#     {#  loop used below is to avoid S memory problems
#      if (verbose) cat("Spawning processes to calculate criteria test for state dimension 1 to ",n)
#      forloop <- function(largeModel, data, warn=TRUE)
#           {if(!is.null(largeModel$G)) z <-largeModel$G[1:forloop.i,,drop=FALSE]
#            else                       z <-NULL
#            z <-SS(F.=largeModel$F[1:forloop.i,1:forloop.i,drop=FALSE],G=z,
#                     H=largeModel$H[  , 1:forloop.i, drop=FALSE],
#                     K=largeModel$K[1:forloop.i,,drop=FALSE])
#            informationTestsCalculations(l(z,data, warn=warn))
#           }
#	assign("balance.forloop", forloop, where=1)
#	assign("forloop.n", n,   where=1 )
#	assign("forloop.values", matrix(NA,n,12),   where=1 )
#	assign("forloop.largeModel", largeModel, where=1)
#	assign("forloop.data", data,   where=1 )
#	assign("forloop.warn", warn,   where=1 )
#	on.exit(remove(c("balance.forloop", "forloop.i", "forloop.n", 
#	   "forloop.values", "forloop.largeModel", 
#	   "forloop.data", "forloop.warn"),where=1))
#	For (forloop.i=1:forloop.n,
#	  forloop.values[forloop.i,]<-balance.forloop(forloop.largeModel, 
#         forloop.data, forloop.warn), sync=TRUE)
#       values <-forloop.values
#      }
    dimnames(values) <- list(NULL,c("port","like","aic","bic", 
          "gvc","rice","fpe","taic","tbic","tgvc","trice","tfpe")) 
    if (verbose) cat("\n")
    opt <-apply(values,2,order)[1,]  # minimum
    if (verbose | is.null(criterion))
      {zz<-criteria.table.nheading()
       options(width=120)
       print(values,digits=4)
       cat("opt     ")
       for (i in 1:length(opt)) cat(opt[i],"    ")
       cat("\n")
       zz<-criteria.table.legend()
      }
    if (is.null(criterion))
      {n <- read.int("Enter the state dimension (enter 0 to stop): ")
       if( n<1) stop("TERMINATED! STATE DIMENSION MUST BE GREATER THAN 0!")
      }
    else { n <- opt[criterion == dimnames(values)[[2]]]  }
    if(m==0) z <-NULL 
      else   z <-largeModel$G[1:n,,drop=FALSE]
    model <- SS(description="nested model a la Mittnik",
          F.=largeModel$F[1:n,1:n,drop=FALSE],G=z,
          H=largeModel$H[,1:n,drop=FALSE],K= largeModel$K[1:n,,drop=FALSE], 
          names=seriesNames(data))
    l(model,data, warn=warn)          
}



############################################################

#     model balancing functions   <<<<<<<<<<

############################################################

balanceMittnik <- function(model, n=NULL){
  if (is.TSestModel(model)) model <- TSmodel(model)
  if  (!is.TSmodel(model)) stop("balanceMittnik expecting a TSmodel.")
  m <- nseriesInput(model)
  newmodel <- SVDbalanceMittnik(markovParms(model), m, n)$model
  newmodel$description <- paste(model$description,"converted to", newmodel$description)
  seriesNames(newmodel)  <-  seriesNames(model)
  newmodel
}


SVDbalanceMittnik <- function(M, m, n=NULL)
{#                   # Form k block Hankel Matrix from M.

   read.int <- function(prmt)
     {err <-TRUE
      while (err)
       {cat(prmt)
 	n <- as.integer(readline()) # crude. this truncates reals
 	if (is.na(n)) cat("value must be an integer\n")
 	else err <- FALSE
       }
      n
     }

   p <- nrow(M)     # dim of endo. series
   r <- m + p       # each sub-matrix is p x r   (= p x (m+p) )
   k<- dim(M)[2] / r
   Hkl <- matrix(0, p*k, r*k) 
   for(i in 1:k)             # Hankel with 0s in SE half  
     Hkl[(1+p*(i-1)):(p*i), 1:(r*(1+k-i))]<- M[,(1+r*(i-1)):(r*k),drop=FALSE]
# note: if last block of M is 0, which is often the case, then filling
#   with zeros, as above, is correct.
#   k <- k %/% 2
#   Hkl <- matrix(0, p*k, r*k) 
#   for (i in 1:k)              # (smaller) completely filled Hankel
#     {for (j in 1:k)       
#        {Hkl[(1+p*(i-1)):(p*i), (1+r*(j-1)):(r*j)] <-M[ ,(1+r*(i+j-2)):(r*(i+j-1))] }}

   svd.of.Hkl <- La.svd(Hkl)
   shifted.Hkl <- Hkl[,(r+1):dim(Hkl)[2]]
   shifted.Hkl <- cbind(shifted.Hkl,matrix(0,p*k,r))
   rtQ <- diag(sqrt(svd.of.Hkl$d))
   rtQ.inv <- diag(1/sqrt(svd.of.Hkl$d))   
   o.mtr <- svd.of.Hkl$u %*% rtQ 
   o.inv <- t(svd.of.Hkl$u %*% rtQ.inv )
   c.mtr <- rtQ %*% svd.of.Hkl$vt
   c.inv <- t(rtQ.inv %*% svd.of.Hkl$vt)
   svd.crit <- svd.of.Hkl$d
   c.mtr[svd.crit==0,] <-0     # this gets rid of NAs but the resulting
   c.inv[,svd.crit==0] <-0     #  model may not be minimal
   o.inv[svd.crit==0,] <-0     # eg. F may have zero rows and columns
   o.mtr[,svd.crit==0] <-0     #   (specifically if the model included a trend)
   if (is.null(n))
     {svd.crit <- svd.criteria(svd.of.Hkl$d)
      n <- read.int("Enter the number of singular values to use for balanced model: ")
     }
   H  <- o.mtr[1:p,1:n,drop=FALSE]
   FF <- o.inv[1:n,,drop=FALSE] %*% shifted.Hkl%*%c.inv[,1:n,drop=FALSE]
   if (m == 0) G <- NULL             # no exog.
   else        G <- c.mtr[1:n,1:m,drop=FALSE]  # exog.
   K  <- c.mtr[1:n,(m+1):(m+p),drop=FALSE]
   FF <- FF + K%*%H  #converts from system with lagged outputs as 
                    #   inputs see Mittnik p1190.
  #browser()
  #  good checks here are (with n=length(svd.of.Hkl$d)):
  #   0 == max(abs(shifted.Hkl - o.mtr[,1:n] %*% (FF-K%*%H) %*% c.mtr[1:n,]))
  #   0 != max(abs(shifted.Hkl))
  model<-SS(description="nested-balanced model a la Mittnik",
              F.=FF,G=G,H=H,K= K)
  list(crit=svd.crit,model=model)
}


MittnikReducedModels <- function(largeModel)
{# Return a list of models with all smaller state dimesions.
  largeModel <- TSmodel(toSS(largeModel))
  largeModel <- balanceMittnik(largeModel, n=dim(largeModel$F)[1])
  r <- vector("list", dim(largeModel$F)[1])
  for (j in 1:length(r))
    r[[j]] <- SS(F.=largeModel$F[1:j,1:j,drop=FALSE],
                G=if(is.null(largeModel$G)) NULL else largeModel$G[1:j,,drop=FALSE],
                H=largeModel$H[  , 1:j, drop=FALSE],   K= largeModel$K[1:j,,drop=FALSE])
  r
}


############################################################################
#
#       experimental estimation techniques    <<<<<<<<<<
#
############################################################################

estBlackBox2 <- function(data, estimation="estVARXls", 
          lag.weight=.9, 
          reduction="MittnikReduction", 
          criterion="taic", 
          trend=FALSE, 
          subtract.means=FALSE,  re.add.means=TRUE, 
          standardize=FALSE, verbose=TRUE, max.lag=12)
{if ((estimation!="estVARXls") && (trend) )
     {cat("Trend estimation only support with estVARXls.\n")
      cat("Proceeding using estVARXls.\n")
      estimation<-"estVARXls"
     }

 if(estimation=="estVARXls")
     model <- estVARXls(data,trend=trend, subtract.means=subtract.means, 
                         re.add.means=re.add.means, max.lag=max.lag, 
                         standardize=standardize, lag.weight=lag.weight)
 # else if(estimation=="estVARXar")
 #     model <-estVARXar(data, subtract.means=subtract.means, max.lag=max.lag)
 else
    stop("Only estVARXls estimation technique is supported to date.\n")
 if (verbose) cat("First VAR model,              lags= ", 
          dim(model$model$A)[1]-1, ", -log likelihood = ", model$estimates$like[1], "\n")
 model <- l(toSS(model),model$data) # data is standardized if standardize=T in estimation
 n <- dim(model$model$F)[1]
 if (verbose) cat("Equivalent    state space model, n= ", 
                  n, ", -log likelihood = ", model$estimates$like[1], "\n")
 if (1 < n)
   {model <- do.call(reduction,
                     list(model,criterion=criterion, verbose=verbose))
   #model <- eval(call(reduction,model,criterion=criterion, verbose=verbose))
    if (verbose) cat("Final reduced state space model, n= ",
              dim(model$model$F)[1], ", -log likelihood = ", model$estimates$like[1], "\n")
   }
  if (verbose &&  dev.cur() != 1 ) checkResiduals(model)
 model
}


bestTSestModel <- function(models, sample.start=10, sample.end=NULL,
    criterion="aic", verbose=TRUE){
  values <- NULL
  for (lst in models ) 
    {z <- informationTestsCalculations(lst, sample.start=sample.start, 
                  sample.end=sample.end)
     values <-rbind(values,z)
    }
  if (verbose)
     {cat("Criterion value for all models based on data starting in period: ",  
           sample.start, "\n")
      cat(values[,criterion], "\n")
     }
#Rbug rbind above looses dimnames
  dimnames(values) <- dimnames(z)
  opt <-order(values[,criterion])[1]  # minimum
  invisible(models[[opt]])
}


estBlackBox3 <- function(data, estimation="estVARXls", 
       lag.weight=1.0, 
       reduction="MittnikReduction", 
       criterion="aic", 
       trend=FALSE, subtract.means=FALSE,  re.add.means=TRUE, 
       standardize=FALSE, verbose=TRUE, max.lag=12, sample.start=10) {
    if ((estimation!="estVARXls") && (trend) )
     {cat("Trend estimation only support with estVARXls.\n")
      cat("Proceeding using estVARXls.\n")
      estimation<-"estVARXls"
     }
 models <- vector("list", max.lag)
 for (i in 1:max.lag)
   {if(estimation=="estVARXls")
      models[[i]] <- estVARXls(data,trend=trend, 
                          subtract.means=subtract.means,
                          re.add.means=re.add.means, max.lag=i, 
                          standardize=standardize, lag.weight=lag.weight)
    else
      stop("Only estVARXls estimation technique is supported to date.\n")
   }
 model <- bestTSestModel(models, criterion=criterion, sample.start=sample.start, verbose=verbose)
 if (verbose) cat("Selected VAR model,              lags= ", 
          dim(model$model$A)[1]-1, ", -log likelihood = ", model$estimates$like[1], "\n")
 model <- l(toSS(model),model$data) # data is standardized if standardize=T in estimation
 n <- dim(model$model$F)[1]
 if (verbose) cat("Equivalent    state space model,    n= ", 
                  n, ", -log likelihood = ", model$estimates$like[1], "\n")
 if (1 < n)
   {model <- do.call(reduction,
                      list(model,criterion=criterion, verbose=verbose))
   #model <- eval(call(reduction,model,criterion=criterion, verbose=verbose))# , sample.start=sample.start))
    if (verbose) cat("Final reduced state space model, n= ",
              dim(model$model$F)[1], ", -log likelihood = ", model$estimates$like[1], "\n")
   }
  if (verbose &&  dev.cur() != 1 ) checkResiduals(model)
 model
}


bft <- function(data, ...) estBlackBox4(data, ...)

estBlackBox4 <- function(data, estimation="estVARXls", 
                lag.weight=1.0,  variable.weights=1, 
                reduction="MittnikReduction", 
                criterion="taic", 
                trend=FALSE, subtract.means=FALSE,  re.add.means=TRUE, 
                standardize=FALSE, verbose=TRUE, max.lag=12, 
		sample.start=10, warn=TRUE)
{if ((estimation!="estVARXls") && (trend) )
     {cat("Trend estimation only support with estVARXls.\n")
      cat("Proceeding using estVARXls.\n")
      estimation<-"estVARXls"
     }
 models <- vector("list", max.lag)
 for (i in 1:max.lag)
   {if(estimation=="estVARXls")
      {model<- estVARXls(data,trend=trend, 
                          subtract.means=subtract.means,
                          re.add.means=re.add.means, max.lag=i, 
                          standardize=standardize, lag.weight=lag.weight,
                          warn=warn)
      }
    else if(estimation=="estWtVariables")
      {model<- estWtVariables(data, variable.weights,
                        estimation.args=list(trend=trend, 
                          subtract.means=subtract.means,
                          re.add.means=re.add.means, max.lag=i, 
                          standardize=standardize, lag.weight=lag.weight,
                          warn=warn) )
      }
    else
      stop("Estimation technique not yet is supported.\n")
    if (verbose) cat("Estimated  VAR   model       -log likelihood = ", 
        model$estimates$like[1],", lags= ",  dim(model$model$A)[1]-1,"\n")
    model <- l(toSS(model),model$data,warn=warn) # data is standardized if standardize=T in estimation
    n <- dim(model$model$F)[1]
    if (verbose) cat("Equivalent state space model -log likelihood = ", 
               model$estimates$like[1], ",   n = ", n, "\n")
    if (1 < n)
      {model <- do.call(reduction,
                    list(model,criterion=criterion, verbose=verbose, warn=warn))
      #model <- eval(call(reduction,model,criterion=criterion, 
      #              verbose=verbose, warn=warn))# , sample.start=sample.start))
       if (verbose) cat("Final reduced state space model, n= ",
           dim(model$model$F)[1], ", -log likelihood = ", 
                         model$estimates$like[1], "\n")
      }
     models[[i]] <- model
   }
 model <- bestTSestModel(models, criterion=criterion, sample.start=sample.start, verbose=verbose)
 if (verbose &&  dev.cur() != 1 ) checkResiduals(model)
 model
}

#  z<-estBlackBox4(eg.data,  max.lag=3 ) 



############################################################

#     statistical test functions  <<<<<<<<<<

############################################################

Portmanteau <- function(res){
  # Portmanteau statistic for residual
  # before R 1.9.0 required ts not stats
  ac <- stats::acf(as.ts(res),type="covariance", plot=FALSE)$acf
  p <- dim(ac)[1]
#  a0 <- solve(ac[1,,])  the following is more robust than solve for
#              degenerate densities
  v <- La.svd(ac[1,,])
#  if(1 == length(v$d)) a0 <- t(v$vt) %*%     (1/v$d) %*% t(v$u)
#	           else a0 <- t(v$vt) %*% diag(1/v$d) %*% t(v$u)	
#  following is equivalent
  a0 <-  t(v$vt) %*% sweep(t(v$u),1,1/v$d, "*") 
  P <-0
  for (i in 2:p) 
    {P <- P + sum(diag(t(ac[i,,]) %*% a0 %*% ac[i,,] %*% a0))
    }
  dim(res)[1]*P
}


checkResiduals <- function(obj, ...)  UseMethod("checkResiduals")

checkResiduals.TSestModel <- function(obj, ...){
   invisible(checkResiduals(obj$estimates$pred - outputData(obj), ...))}

checkResiduals.TSdata <- function(obj, ...){
	invisible(checkResiduals(outputData(obj), ...)) }

checkResiduals.default <- function(obj, ac=TRUE, pac=TRUE,
     select=seq(nseries(obj)), drop=NULL, plot.=TRUE, 
     graphs.per.page=5, verbose=FALSE, ...)
{#  (... further arguments, currently disregarded)
 #  select can indicate column indices of  residuals to use.
 # If drop is not null it should be a vector of the row indices of residuals
 #  which are not to be used. Typically this can be used to get rid
 #  of bad initial conditions (eg. drop=seq(10) ) or outliers.
  resid <- selectSeries(obj, series=select)
  if (!is.null(drop))  resid <- resid[-drop,, drop=FALSE]
  mn<-colMeans(resid)
  acr <-NULL
  pacr<-NULL
  cov <- var(resid)
  if (verbose) 
    {cat("residual means: \n") ; print(mn); cat("\n")
     cat("residual cov  : \n") ; print(cov) ; cat("\n")
    }
  p <- nseries(resid)
  resid0 <- sweep(resid, 2, mn, FUN="-")
#  resid0 <- resid - t(array(colMeans(resid),rev(dim(resid)))) # mean 0
  cusum <- apply(resid0,2,cumsum)/ t(array(diag(var(resid0)),rev(dim(resid0))))
  # before R 1.9.0 required ts not stats
  if(plot. &&  dev.cur() != 1 ) 
    {graphs.per.page <- min(p, graphs.per.page)
     names <- seriesNames(resid)
     old.par <-par(mfcol = c(3, graphs.per.page),
           mar = c(5.1, 4.1,3.1, 0.1), no.readonly=TRUE ) #c(5,4.1,5,0.1) c(2.1, 4.1,3.1, 0.1)
     on.exit(par(old.par))
     for (i in 1:p) 
          {tfOnePlot(resid[,i], xlab=names[i])
           if (i %% graphs.per.page == ceiling(graphs.per.page/2)) 
                 title(main ="Residuals ")
           if (i %% graphs.per.page == 0) 
                 title(main = paste("       page ", floor(i / graphs.per.page)))
           tfOnePlot(cusum[,i], xlab=names[i])
           if (i %% graphs.per.page == ceiling(graphs.per.page/2)) 
                 title(main = "Cusum")
           if (exists("density"))
                       rd <- density(resid[,i],       bw=var(resid[,i])^0.5)
           else if (exists("ksmooth") & ! is.R()) 
                       rd<-ksmooth(resid[,i],bandwidth=var(resid[,i])^0.5)
           else
        stop("Neither ksmooth nor density are available to calculate the plot.")
           plot(rd,type="l",xlab=names[i],ylab="")
           if (i %% graphs.per.page == ceiling(graphs.per.page/2))
              title(main = "kernel estimate of residual distribution")
          }
     if (ac)
       {#par(mfrow = c(1, 1), mar = c(2.1, 4.1,3.1, 0.1), oma=c(0,0,5,0) )
        #acr  <-acf(as.ts(resid), plot=TRUE)$acf
        #mtext("Autocorrelations", side=3, outer=TRUE, cex=1.5)
        acr  <-acf(as.ts(resid), plot=TRUE, mar=c(3,3,2,0.8), oma=c(1,1.2,4,1))$acf
        mtext("Autocorrelations", side=3, outer=TRUE, cex=1.5, line=3)
       }
     if (pac)
       {pacr <- acf(as.ts(resid), plot=TRUE, type= "partial",
	           mar=c(3,3,2,0.8), oma=c(1,1.2,4,1))$acf
        mtext("Partial Autocorrelations", side=3, outer=TRUE, cex=1.5, line=3)
       }
    }
  else 
    {if (ac)  acr  <- acf(as.ts(resid), plot=FALSE)$acf
     if (pac)	  pacr <- acf(as.ts(resid), plot=FALSE, type= "partial")$acf
    }
  if (ac & verbose)
    {cat("residual auto-correlations:\n")
     cat("      lag:   ")
     for (i in 1: dim(acr)[1])  cat(i,"        ")  
     cat("\n") 
     for (i in 1:p) { cat(i,": "); cat(acr[,i,i]); cat("\n")  }
    }
  if (pac & verbose)
    {cat("partial auto-correlations:\n")
     cat("      lag:   ")
     for (i in 1: dim(pacr)[1])  cat(i,"        ")  
     cat("\n") 
     for (i in 1:p) { cat(i,": "); cat(pacr[,i,i]); cat("\n")  }
    }
#  cat("residual normality tests:\n")
#  cat("hetros. tests:\n")

  skewness <- colMeans(sweep(resid,2,mn)^3) / diag(cov)^(3/2)
  kurtosis <- colMeans(sweep(resid,2,mn)^4) / diag(cov)^2
  invisible(list(residuals=resid, mean=mn, cov=cov, acf=acr, pacf=pacr, 
                 cusum=cusum, skewness=skewness, kurtosis=kurtosis))
}

informationTests <- function(..., sample.start=1,sample.end=NULL,
		 Print=TRUE, warn=TRUE){
  if (Print) criteria.table.heading()
  values <- NULL
  options(width=100)
  for (lst in list(...) ) 
    {z <-informationTestsCalculations(lst,
    		 sample.start=sample.start,sample.end=sample.end, warn=warn)
     values <-rbind(values, z)
     if (Print)
       {print(c(z),digits=4)
        cat("\n")
       }
    }
  if (Print & (1 <dim(values)[1]) )
    {cat("opt     ")
     opt <-apply(values,2,order)[1,]  # minimum
     for (i in 1:length(opt)) cat(opt[i],"      ")
     cat("\n")
    }
  if (Print) criteria.table.legend()
  invisible(values)
  }


informationTestsCalculations <- function(lst,
      sample.start=1,sample.end=NULL, warn=TRUE){
   resid <- residuals(lst)
    # the following line is just to work around a bug with old style time series
   if (ncol(outputData(lst$data))==1) dim(resid) <- dim(outputData(lst$data))
   if (is.null(sample.end)) sample.end <- nrow(resid)
   resid <- resid[sample.start:sample.end,,drop=FALSE]
   ml <-   residualStats(resid, NULL, warn=warn)$like[1] # neg.log(likelihood)
# previously   ml <-   L(resid)[1]     # neg. log( likelihood ).
   n  <-  length(coef(lst$model))      # No. of parameters.
   # nt is theorical dimension of parameter space    n(m+2p)
   if (is.ARMA(lst$model)) nt <- NA 
   if (is.SS(lst$model))  
      if (is.null(lst$model$G)) nt <- nrow(lst$model$F)*2*nrow(lst$model$H)
      else   nt <- nrow(lst$model$F)*(ncol(lst$model$G)+2*nrow(lst$model$H))
   r  <- nrow(lst$estimates$pred)*ncol(lst$estimates$pred) #No. of residuals.
   port <-Portmanteau(resid)
   aic <-  2*ml + 2*n             # AIC
   ops  <- options(warn=-1)
   on.exit(options(ops))
   bic <-  2*ml + n*log(r)        # BIC
   gvc <-  2*ml - 2*r*log(1-n/r)  # GCV
   rice<-  2*ml - r*log(1-2*n/r)  # RICE
   fpe <-  2*ml + r*(log(1+(n/r))-log(1-(n/r)))  # FPE
   taic <-  2*ml + 2*nt             # AIC
   tbic <-  2*ml + nt*log(r)        # BIC
   tgvc <-  2*ml - 2*r*log(1-nt/r)  # GCV
   trice<-  2*ml - r*log(1-2*nt/r)  # RICE
   tfpe <-  2*ml + r*(log(1+(nt/r))-log(1-(nt/r)))  # FPE
   z<- matrix(c(port,ml,aic,bic,gvc,rice,fpe,taic,tbic,tgvc,trice,tfpe),1,12)
   dimnames(z)<-list(NULL,c("port","like","aic","bic","gvc","rice","fpe","taic",
                   "tbic","tgvc","trice","tfpe"))
   z
}

criteria.table.heading <- function(){
        cat("                    based on no.of parameters     based on theoretical parameter dim.","\n")
        cat("       PORT  -ln(L)  AIC   BIC   GVC  RICE   FPE   AIC   BIC   GVC   RICE   FPE\n")
}

criteria.table.nheading <- function(){
        cat("                    based on no.of parameters     based on theoretical parameter dim.","\n")
        cat("dim.   PORT  -ln(L)  AIC   BIC   GVC  RICE   FPE   AIC   BIC   GVC   RICE   FPE\n")
}

criteria.table.legend <- function(){
        cat("  PORT  - Portmanteau test                     ") 
        cat("  -ln(L)- neg. log likelihood                  \n")
        cat("  AIC   - neg. Akaike Information Criterion    ")  
        cat("  BIC   - neg. Bayes  Information Criterion    \n")
        cat("  GVC   - Generalized Cross Validation         ") 
        cat("  RICE  - Rice Criterion                       \n")
        cat("  FPE   - Final Prediction Error               \n") 
        cat(" WARNING - These calculations do not account for trend parameters in ARMA models.\n") 
}


svd.criteria <- function(sv){
   cat("\n SINGULAR VALUES OF THE HANKEL MATRIX:\n")
   print(sv,digits=3)
   svsqr.vct <- sv^2/sum(sv^2)    
   cat("\n GELFAND & YAGLOM INFORMATION CRITERIA:\n")
   GY <-cumsum(log(1-svsqr.vct))/sum(log(1-svsqr.vct))
   print(GY,digits=3)
   cat("\n RATIO OF SINGULAR VALUES TO MAX SINGULAR VALUE:\n")
   print(sv/sv[1],digits=3)
   cbind(sv,GY,sv/sv[1])
}   

#######################################################################

tfwindow.TSdata <- function(x, tf=NULL, start=tfstart(tf), end=tfend(tf), warn=TRUE)
{# window a TSdata object
  if (0 != nseriesInput(x))  inputData(x) <- 
      tfwindow(inputData(x), start=start, end=end, tf=tf, warn=warn)
  if (0 != nseriesOutput(x)) outputData(x) <- 
      tfwindow(outputData(x), start=start, end=end, tf=tf, warn=warn)
  x
}

window.TSdata <- function(x, start=NULL, end=NULL, tf=NULL, warn=TRUE, ...)
   tfwindow.TSdata(x, start=start, end=end, tf=tf, warn=warn)
#  (... further arguments, currently disregarded)


combine <- function(e1,e2)UseMethod("combine")
combine.default <- function(e1,e2){list(e1,e2)}


combine.TSdata <- function(e1,e2)
{# make a new TSdata object with the two objects
 outputData(e1) <- tbind(outputData(e1),outputData(e2)) 
 if ((0 != (nseriesInput(e1))) & (0 != (nseriesInput(e2))))
       inputData(e1) <- tbind(inputData(e1),inputData(e2)) 
 else  {if (0 != (nseriesInput(e2)))  inputData(e1) <-inputData(e2)  }
 e1
}


trimNA.TSdata <- function(x, startNAs=TRUE, endNAs=TRUE){
 p <- nseriesOutput(x)
 m <- nseriesInput(x)
 if (m==0)
   mat <- trimNA(outputData(x))
 else
   mat <- trimNA(tbind(inputData(x),outputData(x)),startNAs=startNAs,endNAs=endNAs)
 tf <- tframe(mat)
 if (m!=0)
   {sn <- seriesNamesInput(x)
    inputData(x)  <- tframed(mat[,1:m, drop=FALSE], tf) 
    seriesNamesInput(x) <- sn
   }
 sn <- seriesNamesOutput(x)
 outputData(x) <- tframed(mat[,(m+1):(m+p), drop=FALSE], tf)
 seriesNamesOutput(x) <- sn
 x
}



percentChange.TSestModel <- function(obj, base=NULL, lag=1,
   cumulate=FALSE, e=FALSE, ...)
  {#  (... further arguments, currently disregarded)
      TSdata(output=percentChange(obj$estimates$pred))
  }

percentChange.TSdata <- function(obj, base=NULL, lag=1,
   cumulate=FALSE, e=FALSE, ...)
  {#  (... further arguments, currently disregarded)
   if (0 != (nseriesInput(obj)))  inputData(obj)  <- percentChange(inputData(obj))
   if (0 != (nseriesOutput(obj))) outputData(obj) <- percentChange(outputData(obj))
   obj
  }



############################################################

#     model and data scaling functions   <<<<<<<<<<

############################################################

scale.TSdata <- function(x, center=FALSE, scale=NULL) 
{# scale should be a list with two matrices or vectors, named input and output,
 # giving the multiplication factor for inputs and outputs.
 # vectors are treated as diagonal matrices.
 # If input or output are NULL then no transformation is applied.
 # The resulting data has inputs and outputs which are different from
 #  the original in proportion to scale. ie. if S and T are output and input
 #  scaling matrices then 
 #          y'(t) = S y(t) where y' is the new output
 #          u'(t) = S u(t) where u' is the new input
 sc <- scale$input # inputData(scale) causes problems if scale is a vector
 if (!is.null(sc))
   {if (! (is.matrix(sc) | is.vector(sc))) stop("input scale must be a vector or matrix")
    d <- inputData(x)
    tf <- tframe(d)
    names <- seriesNames(d)
    if(is.matrix(sc))      d <- d %*% t(sc)
    else if(1==length(sc)) d <- d * sc 
    else                   d <- d %*% diag(sc)                             
    tframe(d) <- tf
    seriesNames(d) <- names
    inputData(x) <- d 
   }
 sc <- scale$output # outputData(scale) causes problems if scale is a vector
 if (!is.null(sc))
   {if (! (is.matrix(sc) | is.vector(sc))) stop("output scale must be a vector or matrix")
    d <- outputData(x)
    tf <- tframe(d)
    names <- seriesNames(d)
    if(is.matrix(sc))      d <- d %*% t(sc)
    else if(1==length(sc)) d <- d * sc 
    else                   d <- d %*% diag(sc)                             
    tframe(d) <- tf
    seriesNames(d) <- names
    outputData(x) <- d 
   }
 x
}


scale.TSestModel <- function(x, center=FALSE, scale=NULL)
   {scale(TSmodel(x), center=center, scale=scale)}

checkScale <- function(x, scale) UseMethod("checkScale")

checkScale.TSestModel <- function(x, scale){checkScale(TSmodel(x), scale)}

checkScale.TSmodel <- function(x, scale) 
{# This function only checks for some error conditions.
 # inputData(scale) and outputData(scale) cause problems if scale is a vector
 if (!is.null(scale$input))
   {if (is.matrix(scale$input))
       {if (any(La.svd(scale$input)$d == 0))  
        stop("scale input transformations must be non singular.")
       }
    else if(any(scale$input == 0)) stop("input scale elements must be non zero.")
   }
 if (!is.null(scale$output))
   {if (is.matrix(scale$output))
      {if (any(La.svd(scale$output)$d == 0))  
       stop("output scale transformations must be non singular.")
      }
    else if(any(scale$output == 0)) stop("output scale elements must be non zero.")
   }
 invisible(TRUE)
}


scale.innov <- function(x, center=FALSE, scale=NULL)
{if (!checkScale(x, scale)) stop("scaling error.")
 if (!is.null(scale$input)) 
   {sc <- scale$input
    if(is.vector(sc)) sc <- diag(sc, nseriesInput(x))
    x$G <- x$G %*% solve(sc)
   }
 sc <- scale$output
 if(is.vector(sc)) sc <- diag(sc, nseriesOutput(x))          
 x$H <- sc %*% x$H
 x$K <- x$K %*% solve(sc)
 setTSmodelParameters(x)
}

scale.nonInnov <- function(x, center=FALSE, scale=NULL)
{if (!checkScale(x, scale)) stop("scaling error.")
 if (!is.null(scale$input)) 
   {sc <- scale$input
    if(is.vector(sc)) sc <- diag(sc, nseriesInput(x))
    x$G <- x$G %*% solve(sc)
   }
 sc <- scale$output
 if(is.vector(sc)) sc <- diag(sc, nseriesOutput(x))          
 x$H <- sc %*% x$H
 x$R <- sc %*% x$R %*% solve(sc)
 setTSmodelParameters(x)
}

scale.ARMA <- function(x, center=FALSE, scale=NULL)
{if (!checkScale(x, scale)) stop("scaling error.")
 sc <- scale$output
 if(is.vector(sc)) sc <- diag(sc, nseriesOutput(x))   
 x$A <- polyprod(sc, polyprod(x$A, solve(sc)))
 x$B <- polyprod(sc, polyprod(x$B, solve(sc)))
 if (!is.null(x$C)) 
   {sci <- scale$input
    if (!is.null(sci)) 
       {if(is.vector(sci)) sci <- diag(sci, nseriesInput(x))          
        x$C <- polyprod(x$C, solve(sci))
       }
    x$C <- polyprod(sc, x$C)
   }
 if (!is.null(x$TREND))  x$TREND <- t(sc %*% t(x$TREND))
 setTSmodelParameters(x)
}


############################################################

#     Methods which are generic for models and TSdata    <<<<<<<<<<

# generic methods for TS structures (ie with input and output) <<<<<<<<<<
#                   see also tframe.s                          <<<<<<<<<<

############################################################################
# Tobs and Tobs.default are defined in tfame.s

TobsInput <- function(x)UseMethod( "TobsInput")
TobsOutput <- function(x)UseMethod("TobsOutput")

startInput <- function(x)UseMethod("startInput")
startOutput <- function(x)UseMethod("startOutput")

endInput <- function(x)UseMethod("endInput")
endOutput <- function(x)UseMethod("endOutput")

frequencyInput <- function(x)UseMethod("frequencyInput")
frequencyOutput <- function(x)UseMethod("frequencyOutput")




inputData <- function(x,  series=seqN(nseriesInput(x)))UseMethod("inputData")
outputData <- function(x, series=seqN(nseriesOutput(x)))UseMethod("outputData")

inputData.default <- function(x, series=seqN(nseriesInput(x)))
   selectSeries(x$input, series=series)
outputData.default <- function(x, series=seqN(nseriesOutput(x)))
   selectSeries(x$output, series=series)

"inputData<-" <- function(x,  value)  UseMethod("inputData<-")
"outputData<-" <- function(x,  value) UseMethod("outputData<-")

"inputData<-.default" <- function(x,  value) {x$input <- value; x}
"outputData<-.default" <- function(x,  value){x$output <- value; x}

# The logic (revised as of May26, 1998) is that seriesNames can be an attribute
# of any object (like a matrix) and has been moved to tframe.



 seriesNamesInput <- function(x)UseMethod( "seriesNamesInput")
seriesNamesOutput <- function(x)UseMethod("seriesNamesOutput")

 "seriesNamesInput<-" <- function(x, value)UseMethod( "seriesNamesInput<-")
"seriesNamesOutput<-" <- function(x, value)UseMethod("seriesNamesOutput<-")


############################################################################

#    methods for TSmodels  <<<<<<<<<<
# and check elsewhere too

############################################################################

seriesNames.TSmodel <- function(x)
  {list(input=seriesNamesInput(x), output=seriesNamesOutput(x))}

seriesNamesOutput.TSmodel <- function(x)
 {# return output names if available in the object,
  # otherwise return "out" pasted with integers.
  # Note, there is no TSdata in this case, so tframe is not involved
  if (!is.null(attr(x, "seriesNamesOutput")))
                          return(attr(x, "seriesNamesOutput"))
  else if (0 != nseriesOutput(x)) 
                   return(paste("out", seq(nseriesOutput(x)), sep=""))  
  else return(NULL)
 }

seriesNamesInput.TSmodel <- function(x)
 {# return input names if available in the object,
  # otherwise return "in" pasted with integers.
  if (!is.null(attr(x, "seriesNamesInput")))
                          return(attr(x, "seriesNamesInput"))
  else if (0 != nseriesInput(x)) 
                   return(paste("in", seq(nseriesInput(x)), sep=""))  
  else return(NULL)
 }

"seriesNames<-.TSmodel" <- function(x, value)
   { seriesNamesInput(x) <-  value$input
    seriesNamesOutput(x) <-  value$output
    x
   }

"seriesNamesInput<-.TSmodel" <- function(x,  value)
   {if (!is.null(value) && length(value) != nseriesInput(x))
       stop("model input dimension and number of names do not match.")
    attr(x, "seriesNamesInput")  <- value
    x
   }

"seriesNamesOutput<-.TSmodel" <- function(x,  value)
   {if (!is.null(value) && length(value) != nseriesOutput(x))
       stop("model output dimension and number of names do not match.")
    attr(x, "seriesNamesOutput")  <- value
    x
   }


############################################################################

#    methods for TSestModels  <<<<<<<<<<
# and check elsewhere too   (esp. for start, end and frequency)

############################################################################
start.TSestModel <- function(x, ...)tfstart(x$data)
startInput.TSestModel <- function(x)startInput(x$data)
startOutput.TSestModel <- function(x)startOutput(x$data)

end.TSestModel <- function(x, ...)tfend(x$data)
endInput.TSestModel <- function(x)endInput(x$data)
endOutput.TSestModel <- function(x)endOutput(x$data)

frequency.TSestModel <- function(x, ...)tffrequency(x$data)
frequencyInput.TSestModel <- function(x)frequencyInput(x$data)
frequencyOutput.TSestModel <- function(x)frequencyOutput(x$data)

Tobs.TSestModel <- function(x)Tobs(x$data)
TobsInput.TSestModel <- function(x)TobsInput(x$data)
TobsOutput.TSestModel <- function(x)TobsOutput(x$data)

inputData.TSestModel <- function(x, series=seqN(nseriesInput(x)))
    inputData(x$data, series=series)
outputData.TSestModel <- function(x, series=seqN(nseriesOutput(x)))
    outputData(x$data, series=series)

nseriesInput.TSestModel <- function(x)  nseriesInput(x$data)
nseriesOutput.TSestModel <- function(x) nseriesOutput(x$data)

seriesNamesInput.TSestModel <- function(x)
 {m <- seriesNamesInput(x$model)
  d <- seriesNamesInput(x$data)
  if (!all(m==d)) 
    warning("data and model names do not correspond. Model names returned.")
  m
 }

seriesNamesOutput.TSestModel <- function(x)
 {m <- seriesNamesOutput(x$model)
  d <- seriesNamesOutput(x$data)
  if (!all(m==d)) 
    warning("data and model names do not correspond. Model names returned.")
  m
 }

seriesNames.TSestModel <- function(x)
  {list(input=seriesNamesInput(x), output=seriesNamesOutput(x))}

"seriesNames<-.TSestModel" <- function(x, value)
   { seriesNamesInput(x) <- value$input
    seriesNamesOutput(x) <- value$output
    x
   }


"seriesNamesInput<-.TSestModel" <- function(x, value)
   {seriesNamesInput(x$model) <- value;
    seriesNamesInput(x$data ) <- value;
    x
   }
"seriesNamesOutput<-.TSestModel" <- function(x, value) 
   {seriesNamesOutput(x$model) <- value;
    seriesNamesOutput(x$data ) <- value;
    x
   }


############################################################################

#    methods for TSdata  <<<<<<<<<<
#   (there are methods in other packages too)

############################################################################

is.TSdata <- function(obj) { inherits(obj, "TSdata")}

print.TSdata <- function(x, ...)
 {print.tframed <- function(x, ...)
     {nm <- seriesNames(x)
      tm <- time(x) 
      if (1 != tffrequency(x))
         tm <- paste( floor(tm), "[", 1+ (tm %% 1) * tffrequency(x), "]", sep="")
      dimnames(x) <- list(tm, nm) 
      attr(x,"tframe") <- NULL
      attr(x,"seriesNames") <- NULL
      class(x) <- if (all(class(x) == "tframed")) NULL else class(x)[class(x) != "tframed"]
      print(x,...)
      invisible(x)
     }
  if(0 != (nseriesInput(x)))
     {cat("input data:\n")
      print.tframed(inputData(x), ...)
      if(!is.null(x$input.transformations))
          {cat("input.transformations:\n")
           print(x$input.transformations, ...)
          }
      cat("\n")
     }
  if(0 != (nseriesOutput(x)))
     {cat("output data:\n")
      print.tframed(outputData(x),...)
      if(!is.null(x$output.transformations))
         {cat("output.transformations:\n")
          print(x$output.transformations, ...)
         }
     }
   cat("\n")
   if(!is.null(x$retrieval.date))
      cat("retrieval date: ", x$retrieval.date, "   ")
   if(!is.null(x$source)) 
     {cat("source:\n")
      print(x$source)
     }
  invisible(x)
}

summary.TSdata <- function(object, ...)
  {#  (... further arguments, currently disregarded)
   d <- outputData(object)
   if (!is.tframed(d)) d <- as.ts(d)
   st <- tfstart(d)
   en <- tfend(d)
   fr <- tffrequency(d)
   d <- cbind(d,inputData(object))

   classed(list(  # summary.TSdata constructor
      description=object$description,
      start=st,
      end=en,
      freq=fr,
      sampleT= nrow(outputData(object)),
      p=nseriesOutput(object),
      m=nseriesInput(object),
      ave=colMeans(d),
      max=apply(d,2,max),
      min=apply(d,2,min),
      retrieval.date=object$retrieval.date, 
      source=object$source), "summary.TSdata")
  }

print.summary.TSdata <- function(x, digits=options()$digits, ...)
{#  (... further arguments, currently disregarded)
   if (!is.null(x$description)) cat(x$description,"\n")
   cat("start: ", x$start, " end: ",  x$end," Frequency: ", x$freq,"\n")
   cat("sample length = ",x$sampleT,"\n")
   cat("number of outputs=",x$p, "   number of inputs=",x$m, "\n")
   cat("average :\n")
   print(x$ave)
   cat("max    :\n")
   print(x$max)
   cat("min    :\n")
   print(x$min)
   cat("\n")
   if(!is.null(x$retrieval.date)) 
      cat("retrieval date: ", x$retrieval.date, "   ")
   if(!is.null(x$source)) {cat("source:\n"); print(x$source)  }
   invisible(x)
}


tfplot.TSdata <- function(x, ..., 
        tf=NULL, start=tfstart(tf), end=tfend(tf),
        select.inputs  = seq(length=nseriesInput(x)),
        select.outputs = seq(length=nseriesOutput(x)),
        Title=NULL, xlab=NULL, ylab=NULL, 
	graphs.per.page=5, mar=par()$mar, reset.screen=TRUE)
{# plot input and output data.
 # ... is a list of objects of class TSdata (with similar input
 #       and output dimensions.
 #previously mar = if(is.R()) c(3.1,6.1,1.1,2.1) else c(5.1,6.1,4.1,2.1)
 # Note that using ... like this means it cannot be used to pass additional
 #   arguments to plot, so unfortunately all necessary plot arguments must be 
 #   explicit in the arguments to tfplot.TSdata
 # start is the starting point (date)  and end the ending point for
 # plotting. If not specified the whole sample is plotted.
# output graphics can be paused between pages by setting par(ask=T).
  #x <- freeze(x)
  Ngraphs <- length(select.outputs) + length(select.inputs)
  if(reset.screen)
    {Ngraphs <- min(Ngraphs, graphs.per.page)
     old.par <- par(mfcol = c(Ngraphs, 1), mar=mar, no.readonly=TRUE) # previously c(5.1,6.1,4.1,2.1)) 
     on.exit(par(old.par))
    }
  #if (is.null(Title)) Title <- ""

  if (!is.null(ylab))
       names <- rep(ylab, nseriesInput(x) + nseriesOutput(x))
  else names <-  c(seriesNamesInput(x),  seriesNamesOutput(x))

  if (1 == length(xlab)) xlab <- rep(xlab, length(names))

#  if (0 != length(select.inputs)) 
    {for (i in select.inputs) 
      {j <- 0
       z <- matrix(NA, TobsInput(x), length(append(list(x),list(...))))
       if(mode(i)=="character") i <- match(i, seriesNamesInput(x))
       for (d in append(list(x),list(...)) ) 
         {if (!is.TSdata(d))
            stop("Expecting TSdata objects. Do not truncate argument names as that can cause a problem here.")
          #d <- freeze(d)
          j <- j + 1
          z[, j] <- inputData(d, series = i)
         }
       tframe(z) <-tframe(inputData(x))
       tfOnePlot(z, xlab=xlab[i], ylab=names[i], start=start, end=end) 
       if(!is.null(Title) && (i==1) && (is.null(options()$PlotTitles)
                || options()$PlotTitles)) title(main = Title)
    } }
  for (i in select.outputs) 
    {j <-0
     if(mode(i)=="character") i <- match(i, seriesNamesOutput(x))
     z <- matrix(NA,TobsOutput(x),length(append(list(x),list(...))))
     #z <-NULL 
     for (d in append(list(x),list(...)) ) 
       {if (!is.TSdata(d))
            stop("Expecting TSdata objects. Do not truncate argument names as that can cause a problem here.")
        #d <- freeze(d)
        j <- j+1
        z[,j]<-outputData(d, series=i) 
       }
     tframe(z) <-tframe(outputData(x))
     tfOnePlot(z, xlab=xlab[i], ylab=names[nseriesInput(x) + i],start=start, end=end)
     if(!is.null(Title) && (0 == nseriesInput(x)) && (i==1) && 
      (is.null(options()$PlotTitles) || options()$PlotTitles)) title(main = Title)
    }
  invisible()
}
 
inputData.TSdata <- function(x, series=seqN(nseriesInput(x)))
  {if (is.null(x$input)) NULL  else selectSeries(x$input, series=series) }

outputData.TSdata <- function(x, series=seqN(nseriesOutput(x)))
  {if (is.null(x$output)) NULL  else selectSeries(x$output, series=series) }


"inputData<-.TSdata" <- function(x, value) 
   {cls <- class(x); x$input <- value; class(x) <- cls
    x
   }


"outputData<-.TSdata" <- function(x, value)
   {cls <- class(x); x$output <-value; class(x) <- cls
    x
   }

# Note: series names changed to an attribute of input and output for 
#     data but not for models!!!!

seriesNames.TSdata <- function(x)
 {list(input=seriesNamesInput(x), output=seriesNamesOutput(x))}

 seriesNamesInput.TSdata <- function(x) {seriesNames( inputData(x))}
seriesNamesOutput.TSdata <- function(x) {seriesNames(outputData(x))}

"seriesNames<-.TSdata" <- function(x, value){ 
    seriesNamesInput(x)  <-  value$input
    seriesNamesOutput(x) <-  value$output
    x
   }

"seriesNamesInput<-.TSdata"  <- function(x,  value)
   {seriesNames(inputData(x))  <- value; x }

"seriesNamesOutput<-.TSdata" <- function(x,  value) 
   {seriesNames(outputData(x)) <- value; x }

nseriesInput.TSdata <- function(x)
   {if (is.null(x$input)) 0 else nseries(x$input)}

nseriesOutput.TSdata <- function(x)
   {if (is.null(x$output)) 0 else nseries(x$output)}

start.TSdata <- function(x, ...)
{#  (... further arguments, currently disregarded)
 i  <-  startInput(x)
 o  <- startOutput(x)
 if (((!is.null(o)) & (!is.null(i))) && all(i==o)) return(o)
 else return(c(i,o))
}

startInput.TSdata <- function(x)
 {if (is.null(x$input))  return(NULL)
  else return(tfstart(x$input))
 }

startOutput.TSdata <- function(x)
 {if (is.null(x$output))  return(NULL)
  else return(tfstart(x$output))
 }

end.TSdata <- function(x, ...)
{#  (... further arguments, currently disregarded)
 i  <-  endInput(x)
 o  <- endOutput(x)
 if (((!is.null(o)) & (!is.null(i))) && all(i==o)) return(o)
 return(c(i,o))
}

endInput.TSdata <- function(x)
 {if (is.null(x$input))  return(NULL)
  else return(tfend(x$input))
 }

endOutput.TSdata <- function(x)
 {if (is.null(x$output))  return(NULL)
  else return(tfend(x$output))
 }

frequency.TSdata <- function(x, ...)
{#  (... further arguments, currently disregarded)
 i  <-  frequencyInput(x)
 o  <- frequencyOutput(x)
 if (((!is.null(o)) & (!is.null(i))) && all(i==o)) return(o)
 return(c(i,o))
}

frequencyInput.TSdata <- function(x)
 {if (is.null(x$input))  return(NULL)
  else return(tffrequency(x$input))
 }

frequencyOutput.TSdata <- function(x)
 {if (is.null(x$output))  return(NULL)
  else return(tffrequency(x$output))
 }


Tobs.TSdata <- function(x, ...) Tobs(outputData(x))
#  (... further arguments, currently disregarded)
TobsOutput.TSdata <- function(x)  Tobs(outputData(x))[1]
TobsInput.TSdata  <- function(x)  Tobs( inputData(x))[1]

tbind.TSdata <- function(x, d2, ..., pad.start=TRUE, pad.end=TRUE, warn=TRUE)
 {if( ! (is.TSdata(x) & is.TSdata(d2)))
     stop("tbind requires arguments to be of a similar type (ie. TSdata).")
  if ((0 != nseriesInput(x)) || (0 != nseriesInput(d2)) )
    inputData(x) <- tbind(inputData(x),inputData(d2),
                 ..., pad.start=pad.start, pad.end=pad.end, warn=warn)
  if ((0 != nseriesOutput(x)) || (0 != nseriesOutput(d2)) )
    outputData(x) <- tbind(outputData(x),outputData(d2),
                 ..., pad.start=pad.start, pad.end=pad.end, warn=warn)
  x
 }



testEqual.TSdata <- function(obj1, obj2, fuzz=1e-16)
  {r <- TRUE
   if (r & (!is.null(obj1$input)))
     {if(is.null(obj2$input)) r <- FALSE
      else  r <-testEqual(obj1$input, obj2$input, fuzz=fuzz)
     }
   if (r & (!is.null(obj1$output)))
     {if(is.null(obj2$output)) r <- FALSE
      else r <-testEqual(obj1$output, obj2$output, fuzz=fuzz)
     }
   r
  }




TSdata <- function (data=NULL, ...) UseMethod("TSdata") 
 

#TSdata.default <- function(data=NULL, input=NULL, output=NULL, ...)  
#{if (is.null(data) && (!is.null(input) | !is.null(output) ))
#    {if(!is.null(input) && is.vector(input)) input <- 
#	   tframed(matrix(input, length(input),1), tf=tframe(input),
#	           names=seriesNames(input))
#     if(!is.null(output) && is.vector(output)) output <- 
#	   tframed(matrix(output, length(output),1), tf=tframe(output),
#	           names=seriesNames(output))
#     data <- classed(list(input=input, output=output), "TSdata") # constructor
#  }else 
#     data <- classed(data, "TSdata")   # constructor keeps other list elements
#  
# if(!is.list(data)) stop("TSdata.default could not construct a TSdata format.")	    
# if   ( 0 == nseriesInput(data)   &    0 == nseriesOutput(data) 
#    |  (0 != nseriesOutput(data) && !is.matrix(outputData(data)))
#    |  (0 !=  nseriesInput(data) && !is.matrix( inputData(data))) )
#      stop("TSdata.default could not construct a TSdata format.")
# data
#}

TSdata.default <- function(data=NULL, input=NULL, output=NULL, ...)  
{if (is.null(data) & is.null(input)  & is.null(output) )
      stop("TSdata.default needs data, input, or output specified.")
 if (is.null(data))
    {if(!is.null(input)) input <- 
	   tframed(as.matrix(input), tf=tframe(input), names=seriesNames(input))
     if(!is.null(output)) output <- 
	   tframed(as.matrix(output), tf=tframe(output),names=seriesNames(output))
     data <- classed(list(input=input, output=output), "TSdata") # constructor
  }else 
     data <- classed(data, "TSdata")   # constructor keeps other list elements
  
 if(!is.list(data)) stop("TSdata.default could not construct a TSdata format.")	    
 if   ( 0 == nseriesInput(data)   &    0 == nseriesOutput(data) 
    |  (0 != nseriesOutput(data) && !is.matrix(outputData(data)))
    |  (0 !=  nseriesInput(data) && !is.matrix( inputData(data))) )
      stop("TSdata.default could not construct a TSdata format.")
 data
}

TSdata.TSdata <- function(data, ...)  {data} # already TSdata 

TSdata.TSestModel <- function(data, ...)  {data$data} # extract TSdata 

as.TSdata <- function(d) 
 {# Use whatever is actually $input and $output,
  #strip any other class, other parts of the list, and DO NOT use
  # inputData(d) and outputData(d) which may eg reconstitute something
  TSdata(input=d$input, output=d$output)
 }

"tframe<-.TSdata" <- function(x, value) {
 if (0 != nseriesInput(x))  tframe(inputData(x))  <- value
 if (0 != nseriesOutput(x)) tframe(outputData(x)) <- value
 x
}

tframed.TSdata <- function(x, tf=NULL, names=NULL, ...)  
{# switch to tframe representation
 if(0 != (nseriesOutput(x)))
       outputData(x) <-tframed(outputData(x), tf=tf, names=names$output, ...)
 if (0 != (nseriesInput(x)))
        inputData(x) <-tframed(inputData(x),  tf=tf, names=names$input, ...)
 x
}  


############################################################

#  Utility function for time series noise   <<<<<<<<<<

############################################################

makeTSnoise <- function(sampleT,p,lags,noise=NULL, rng=NULL,
                        Cov=NULL, sd=1, noise.model=NULL,
                        noise.baseline=0,
                        tf=NULL, start=NULL,frequency=NULL)
 {# CAUTION: changes here can affect historical comparisons.
  # noise.baseline is added to noise. It should be either a scalar, a matrix of
  #   the same dimension as noise (or noise generated by noise.model), or a
  #   vector of length equal to the dimension of the noise process (which will
  #   be replicated for all Tobs.)
 if(is.null(rng)) rng <- setRNG() # returns setting so don't skip if NULL
 else        {old.rng <- setRNG(rng);  on.exit(setRNG(old.rng))  }

  if ( (!is.null(noise)) & (!is.null(noise.model)) )
    stop("noise and noise.model cannot both be specified.")

  if(!is.null(noise.baseline) && is.matrix(noise.baseline) &&
    (sampleT < dim(noise.baseline)[1]))
      {warning("sampleT (and start date) for noise adjusted to match noise.baseline")
       sampleT <- dim(noise.baseline)[1]
      }

 # Note: noise is added to initial conditions.
 if (!is.null(noise.model))
   {noise <- outputData(simulate(noise.model, sampleT=sampleT+lags))
    noise <- list(w0=noise[1:lags,,drop=FALSE], w=noise[lags+seq(sampleT),,drop=FALSE])
   }
  if (is.null(noise)) {
    w0 <-matrix(NA,lags,p)
    w <- matrix(NA,sampleT,p)
    if (!is.null(Cov)) {
        if(length(Cov) == 1) Cov <- diag(Cov, p)
	W <- t(chol(Cov))
        w <- t(W %*% t(matrix(rnorm((lags+sampleT)*p),(lags+sampleT),p)))
#  above has unfortunate effect that w0 was not from the first p values below
#   would be better, but changes a lot? of tests
#        w <- t(W %*% matrix(rnorm((lags+sampleT)*p),p, (lags+sampleT)))
	w0 <- w[1:lags,]
	w  <- w[-c(1:lags),]
        }
    else {
       if (length(sd)==1) sd <-rep(sd,p)
       for (i in 1:p)
         {w0[,i] <- rnorm(lags,sd=sd[i])
          w[,i]  <- rnorm(sampleT,sd=sd[i])
         }
       }
    noise <- list(w=w, w0=w0)
    }
  else {
     if (is.null(noise$w0) || is.null(noise$w) )
       stop("supplied noise structure is not correct.")
     }
       
  if(!is.null(noise.baseline))
     {if (is.vector(noise.baseline))
        {if(length(noise.baseline)==1) noise$w <- noise$w + noise.baseline
         else if(length(noise.baseline)==1)
           noise$w <-noise$w + t(array(noise.baseline, rev(dim(noise$w))))
         else stop("noise.baseline vector is not correct length.")
        }
      else  noise$w <- noise$w + noise.baseline
     }

  if(!is.null(tf)) tframe(noise$w) <- tf
  else if(!is.null(start))
      {if (is.null(frequency))
         {frequency <- 1
          warning("start set but frequency not specified. Using frequency=1.")
         }
       else noise$w <-tframed(noise$w, list(start=start, frequency=frequency))
       if (is.tframed(noise.baseline) &&
           testEqual(tframe(noise.baseline),tframe(noise$w)))
           {warning("tframe of noise set to tframe of noise.baseline.")
            tframe(noise$w )<-tframe(noise.baseline)
            if(!all(dimnames(noise$w)[[2]] == dimnames(noise.baseline)[[2]]))
              warning("noise names and noise.baseline names do not correspond.")
           }
      }
  append(noise, list(sampleT=sampleT, rng=rng,
     Cov=Cov, sd=sd, noise.model=noise.model,version=as.vector(version)))
 }