R/dse1.R

Defines functions makeTSnoise tframed.TSdata as.TSdata TSdata.TSestModel TSdata.TSdata TSdata.default TSdata testEqual.TSdata tbind.TSdata TobsInput.TSdata TobsOutput.TSdata Tobs.TSdata frequencyOutput.TSdata frequencyInput.TSdata frequency.TSdata endOutput.TSdata endInput.TSdata end.TSdata startOutput.TSdata startInput.TSdata start.TSdata nseriesOutput.TSdata nseriesInput.TSdata seriesNamesOutput.TSdata seriesNamesInput.TSdata seriesNames.TSdata outputData.TSdata inputData.TSdata tfplot.TSdata print.summary.TSdata summary.TSdata print.TSdata is.TSdata seriesNames.TSestModel seriesNamesOutput.TSestModel seriesNamesInput.TSestModel nseriesOutput.TSestModel nseriesInput.TSestModel outputData.TSestModel inputData.TSestModel TobsOutput.TSestModel TobsInput.TSestModel Tobs.TSestModel frequencyOutput.TSestModel frequencyInput.TSestModel frequency.TSestModel endOutput.TSestModel endInput.TSestModel end.TSestModel startOutput.TSestModel startInput.TSestModel start.TSestModel seriesNamesInput.TSmodel seriesNamesOutput.TSmodel seriesNames.TSmodel seriesNamesOutput seriesNamesInput outputData.default inputData.default outputData inputData frequencyOutput frequencyInput endOutput endInput startOutput startInput TobsOutput TobsInput scale.ARMA scale.nonInnov scale.innov checkScale.TSmodel checkScale.TSestModel checkScale scale.TSestModel scale.TSdata percentChange.TSdata percentChange.TSestModel trimNA.TSdata combine.TSdata combine.default combine window.TSdata tfwindow.TSdata svd.criteria criteria.table.legend criteria.table.nheading criteria.table.heading informationTestsCalculations informationTests checkResiduals.default checkResiduals.TSdata checkResiduals.TSestModel checkResiduals Portmanteau estBlackBox4 bft estBlackBox3 bestTSestModel estBlackBox2 MittnikReducedModels SVDbalanceMittnik balanceMittnik MittnikReduction.from.Hankel MittnikReduction estSSMittnik estBlackBox1 estBlackBox estMaxLik.TSmodel estMaxLik.TSdata estMaxLik.TSestModel estMaxLik estWtVariables estSSfromVARX estVARXar estVARXmean.correction fake.TSestModel.missing.data estVARXls state smoother.nonInnov smoother.TSmodel smoother.TSestModel smoother l.SS l.ARMA l.TSestModel l.TSdata l sumSqerror residualStats simulate.ARMA simulate.SS simulate.TSestModel simulate printTestValue DSE.ar setArrays.ARMA setArrays.SS setArrays.TSestModel setArrays setTSmodelParameters.ARMA setTSmodelParameters.SS setTSmodelParameters.TSestModel setTSmodelParameters is.ARMA is.nonInnov.SS is.innov.SS is.SS is.TSestModel is.TSmodel coef.TSestModel coef.TSmodel SS ARMA TSmodel.TSestModel TSmodel.TSmodel TSmodel TSestModel.TSestModel TSestModel checkConsistentDimensions.default checkConsistentDimensions.ARMA checkConsistentDimensions.SS checkConsistentDimensions.TSestModel checkConsistentDimensions.TSdata checkConsistentDimensions nseriesOutput.TSestModel nseriesOutput.ARMA nseriesOutput.SS nseriesOutput.default nseriesOutput nseriesInput.TSestModel nseriesInput.ARMA nseriesInput.SS nseriesInput.default nseriesInput nstates.TSestModel nstates.ARMA nstates.SS nstates companionMatrix characteristicPoly polyvalue polydet polyrootDet polysum polyprod markovParms Riccati toARMA.SS toARMA.ARMA toARMA.TSestModel toARMA toSSChol.TSmodel toSSChol.TSestModel toSSChol fixF toSSOform.TSmodel toSSOform.TSestModel toSSOform toSSinnov fixConstants gmap toSSaugment.ARMA toSSaugment.TSestModel toSSaugment toSSnested.ARMA toSSnested.SS toSSnested.TSestModel toSSnested toSS.ARMA toSS.SS toSS.TSestModel toSS checkBalanceMittnik.ARMA checkBalanceMittnik.SS checkBalanceMittnik.TSestModel checkBalanceMittnik checkBalance.ARMA checkBalance.SS checkBalance.TSestModel checkBalance reachability.ARMA reachability.SS reachability.TSestModel reachability observability.ARMA observability.SS observability.TSestModel observability addPlotRoots plot.roots roots.ARMA roots.SS roots.TSestModel roots stability.ARMA stability.roots stability.TSmodel stability.TSestModel stability McMillanDegree.SS McMillanDegree.ARMA McMillanDegree.TSestModel McMillanDegree testEqual.SS testEqual.ARMA testEqual.TSmodel testEqual.TSestModel tfplot.TSestModel print.summary.ARMA summary.ARMA print.summary.SS summary.SS print.summary.TSestModel summary.TSestModel print.ARMA print.SS print.TSestModel acf.TSestModel acf.TSdata residuals.TSestModel DSEversion onAttach

Documented in acf.TSdata acf.TSestModel addPlotRoots ARMA as.TSdata balanceMittnik bestTSestModel bft characteristicPoly checkBalance checkBalance.ARMA checkBalanceMittnik checkBalanceMittnik.ARMA checkBalanceMittnik.SS checkBalanceMittnik.TSestModel checkBalance.SS checkBalance.TSestModel checkConsistentDimensions checkConsistentDimensions.ARMA checkConsistentDimensions.default checkConsistentDimensions.SS checkConsistentDimensions.TSdata checkConsistentDimensions.TSestModel checkResiduals checkResiduals.default checkResiduals.TSdata checkResiduals.TSestModel checkScale checkScale.TSestModel checkScale.TSmodel coef.TSestModel coef.TSmodel combine combine.default combine.TSdata companionMatrix criteria.table.heading criteria.table.legend criteria.table.nheading DSE.ar DSEversion endInput endInput.TSdata endInput.TSestModel endOutput endOutput.TSdata endOutput.TSestModel end.TSdata end.TSestModel estBlackBox estBlackBox1 estBlackBox2 estBlackBox3 estBlackBox4 estMaxLik estMaxLik.TSdata estMaxLik.TSestModel estMaxLik.TSmodel estSSfromVARX estSSMittnik estVARXar estVARXls estVARXmean.correction estWtVariables fake.TSestModel.missing.data fixConstants fixF frequencyInput frequencyInput.TSdata frequencyInput.TSestModel frequencyOutput frequencyOutput.TSdata frequencyOutput.TSestModel frequency.TSdata frequency.TSestModel gmap informationTests informationTestsCalculations inputData inputData.default inputData.TSdata inputData.TSestModel is.ARMA is.innov.SS is.nonInnov.SS is.SS is.TSdata is.TSestModel is.TSmodel l l.ARMA l.SS l.TSdata l.TSestModel makeTSnoise markovParms McMillanDegree McMillanDegree.ARMA McMillanDegree.SS McMillanDegree.TSestModel MittnikReducedModels MittnikReduction MittnikReduction.from.Hankel nseriesInput nseriesInput.ARMA nseriesInput.default nseriesInput.SS nseriesInput.TSdata nseriesInput.TSestModel nseriesOutput nseriesOutput.ARMA nseriesOutput.default nseriesOutput.SS nseriesOutput.TSdata nseriesOutput.TSestModel nstates nstates.ARMA nstates.SS nstates.TSestModel observability observability.ARMA observability.SS observability.TSestModel outputData outputData.default outputData.TSdata outputData.TSestModel percentChange.TSdata percentChange.TSestModel plot.roots polydet polyprod polyrootDet polysum polyvalue Portmanteau print.ARMA print.SS print.summary.ARMA print.summary.SS print.summary.TSdata print.summary.TSestModel printTestValue print.TSdata print.TSestModel reachability reachability.ARMA reachability.SS reachability.TSestModel residualStats residuals.TSestModel Riccati roots roots.ARMA roots.SS roots.TSestModel scale.ARMA scale.innov scale.nonInnov scale.TSdata scale.TSestModel seriesNamesInput seriesNamesInput.TSdata seriesNamesInput.TSestModel seriesNamesInput.TSmodel seriesNamesOutput seriesNamesOutput.TSdata seriesNamesOutput.TSestModel seriesNamesOutput.TSmodel seriesNames.TSdata seriesNames.TSestModel seriesNames.TSmodel setArrays setArrays.ARMA setArrays.SS setArrays.TSestModel setTSmodelParameters setTSmodelParameters.ARMA setTSmodelParameters.SS setTSmodelParameters.TSestModel simulate simulate.ARMA simulate.SS simulate.TSestModel smoother smoother.nonInnov smoother.TSestModel smoother.TSmodel SS stability stability.ARMA stability.roots stability.TSestModel stability.TSmodel startInput startInput.TSdata startInput.TSestModel startOutput startOutput.TSdata startOutput.TSestModel start.TSdata start.TSestModel state summary.ARMA summary.SS summary.TSdata summary.TSestModel sumSqerror SVDbalanceMittnik svd.criteria tbind.TSdata testEqual.ARMA testEqual.SS testEqual.TSdata testEqual.TSestModel testEqual.TSmodel tfplot.TSdata tfplot.TSestModel tframed.TSdata tfwindow.TSdata toARMA toARMA.ARMA toARMA.SS toARMA.TSestModel TobsInput TobsInput.TSdata TobsInput.TSestModel TobsOutput TobsOutput.TSdata TobsOutput.TSestModel Tobs.TSdata Tobs.TSestModel toSS toSS.ARMA toSSaugment toSSaugment.ARMA toSSaugment.TSestModel toSSChol toSSChol.TSestModel toSSChol.TSmodel toSSinnov toSSnested toSSnested.ARMA toSSnested.SS toSSnested.TSestModel toSSOform toSSOform.TSestModel toSSOform.TSmodel toSS.SS toSS.TSestModel trimNA.TSdata TSdata TSdata.default TSdata.TSdata TSdata.TSestModel TSestModel TSestModel.TSestModel TSmodel TSmodel.TSestModel TSmodel.TSmodel window.TSdata

###
###  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 %*%