R/EvalEst.R

In EvalEst: Dynamic Systems Estimation - Extensions

############################################################################
# Comments here need to be cleaned up to reflect re-org
# Functions in the next group are mainly for evaluating estimation techniques.

# The first group are for generating simulations (ie- generate multiple
#   stochastic simulations of a model using simulate.)
#   These are good for looking at the stochastic properties of a model,
#   but mostly these are intended
#   only for verification purposes since other functions also generate 
#   simulations and it is usually more efficient to regenerate by setting
#   the RNG and seed than it is to save data.
#   The main function in this group is MonteCarloSimulations().
#   The main object class in this group is "simulation"

# The second group are for analysing the convergence of estimators. This
#  is extended to functions of estimates (such as model roots). These
#  functions evaluate a single given estimation method with multiple
#  simulated data sets from a given "true" model.
#  The main function in this group is EstEval().
#  The main object classes in this group are
#    "EstEval"
#     c("rootsEstEval","EstEval")
#     c("TSmodelEstEval","EstEval")
#     c("TSestModelEstEval","EstEval")
#     c("coefEstEval","EstEval")           and
#     c("rootsEstEval","EstEval")


# The third group applies multiple estimation techniques to a given data set.
#  This is primarily a utility for other functions 
#  The main (only) function in this group is estimateModels().
#  It returns an object of class c("estimatedModels")


# The fourth group looks at the forecasts and the covariance of forecasts 
#   for multiple horizons. 
#  The simplest case horizonForecasts() which calculates the forecast for 
#    different horizons at all periods in the sample. This is primarily
#    a utility for calculating the forecast error.
#  The next case is is estimatorsHorizonForecastsWRTdata() which 
#   is an extention of horizonForecasts(). It takes specified data 
#   and estimation techniques and calculates forecasts from the estimated 
#   models.
#  The generic function forecastCov() which considers mulitple  
#   models and calculates the cov of forecasts relative to a given data set.
#   It takes a list of models (+ trend, zero)  and calculates the cov 
#   of predictions. It uses forecastCovSingleModel()
#  The next case, forecastCovEstimatorsWRTdata() uses a list of estimation
#   methods to estimate a list of models and calculate the cov of predictions 
#   relative to one given data set.

#   The next case forecastCovWRTtrue() takes a list of models (+ trend,
#   zero)  and calculates the cov of forecasts relative to data sets 
#   simulated with a true.model.
#   The next case, forecastCovEstimatorsWRTtrue simulates data and 
#   uses a list of estimation methods to estimate a list of models, then
#   calculates the cov of predictions relative to other simulated data set.
#  The main object classes in this group are
#     c("estimatorsHorizonForecastsWRTdata") # ? "horizonForecasts")
#     "horizonForecasts"
#     c("multiModelHorizonForecasts","horizonForecasts")
#     "forecastCov"
#     c("forecastCovWRTdata", "forecastCov")
#     c("forecastCovWRTtrue", "forecastCov")
#     c("forecastCovEstimatorsWRTdata",  "forecastCov")
#     c("forecastCovEstimatorsWRTtrue",  "forecastCov")


# The fifth group are some experimental estimation techniques.

############################################################################
#
#       methods for MonteCarloSimulations  <<<<<<<<<<
#
############################################################################

generateSSmodel <- function(m,n,p, stable=FALSE)
 {#randomly generate an innov state space model. Discard models with largest root
  # greater than 1 (if stable=F) or equal to or greater than 1 if stable=T.
  repeat
    {FF <- matrix(runif(n^2, min=-1,max=1),n,n)
     if (m!=0) G <- matrix(runif(n*m, min=-1,max=1),n,m)
     else G <- NULL
     H <- matrix(runif(n*p, min=-1,max=1),p,n)
     K <- matrix(runif(n*p, min=-1,max=1),n,p)
     model <- SS(F.=FF, G=G,H=H,K=K)
     if (stable) {if (max(Mod(roots(model))) <  1.0) break()}
     else        {if (max(Mod(roots(model))) <= 1.0) break()}
    }
  model
 }


MonteCarloSimulations <- function(model, simulation.args=NULL, 
           replications=100, rng=NULL, quiet=FALSE, ...) UseMethod("MonteCarloSimulations")

MonteCarloSimulations.default <- function (model, simulation.args=NULL, 
 		replications=100, rng=NULL, quiet=FALSE, ...){
        #  (... further arguments, currently disregarded)
 	if (is.null(rng)) rng <- setRNG()
	else {
		old.rng <- setRNG(rng)
		on.exit(setRNG(old.rng))
	}
	arglist <- append(list(model), simulation.args)
	r <- do.call("simulate", arglist)
	if (! is.matrix(r)) stop("simulate(model) must return a matrix.")
        result <- array(NA, c(dim(r), replications))
        tfr <- tframe(r)
        result[, , 1] <- r
        if (1 < replications) 
            for (i in 2:replications)
	      result[, , i] <- do.call("simulate", arglist)
	# default does not work on array seriesNames(result) <- seriesNames(r)
        if (length(seriesNames(r)) != dim(result)[2])
         stop("length of names (",length(seriesNames(r)),
	      ") does not match number of series(",dim(result)[2],").")
        attr(result,"seriesNames") <- seriesNames(r)
	result <- tframed(result, tfr)
	invisible(classed(list(simulations = result, model = model, 
        rng = rng,  simulation.args = simulation.args, 
        description = "data generated by MonteCarloSimulations.default"), 
        c("MonteCarloSimulations")))
}


MonteCarloSimulations.TSestModel <- function(model, simulation.args=NULL, 
           replications=100, rng=NULL, quiet=FALSE, ...)
  {if (is.null(simulation.args$sd) & is.null(simulation.args$Cov)) 
     simulation.args$Cov <- model$estimates$cov
   if (is.null(simulation.args$input)) simulation.args$input <- inputData(model)
   MonteCarloSimulations(TSmodel(model), simulation.args=simulation.args, 
           replications=replications, rng=rng, quiet=quiet, ...)
  }

MonteCarloSimulations.EstEval <- function(model, simulation.args=NULL,
            replications=100, rng=getRNG(model),  quiet=FALSE, ...)
     MonteCarloSimulations(TSmodel(model), simulation.args=simulation.args,
         replications=replications, rng=rng,  quiet=quiet, ...)
# this looks like a candidate for NextMethod

MonteCarloSimulations.MonteCarloSimulations <- function(model, simulation.args=NULL,
            replications=100, rng=getRNG(model),  quiet=FALSE, ...)
     MonteCarloSimulations(TSmodel(model), simulation.args=simulation.args,
         replications=replications, rng=rng,  quiet=quiet, ...)
# this looks like a candidate for NextMethod


MonteCarloSimulations.TSmodel <- function(model, simulation.args=NULL,
          replications=100, rng=NULL, quiet=FALSE, ...)
#	  Spawn=if (exists(".SPAWN")) .SPAWN else FALSE, ...)
{ 
 if(is.null(rng)) rng <- setRNG() # returns setting so don't skip if NULL
 else        {old.rng <- setRNG(rng);  on.exit(setRNG(old.rng))  }
 
 arglist <- append(list(model), simulation.args)
# if (Spawn)
#  {if (!quiet)cat("Spawning processes to calculate ", replications, " replications.\n")
#   assign("sim.forloop.n", replications, where = 1)
#   assign("sim.forloop.result", list(NULL), where = 1)
# #  assign("sim.forloop.model", model, where = 1)
#   assign("sim.forloop.arglist", arglist, where = 1)
#   on.exit(remove(c("sim.forloop.i", "sim.forloop.n", "sim.forloop.result",
#       "sim.forloop.arglist"),where = 1))
#   For(sim.forloop.i = 1:sim.forloop.n, sim.forloop.result[[sim.forloop.i]] <- 
#       do.call("simulate",  sim.forloop.arglist),
#       first=options(warn=-1), sync = TRUE)
#   result <- array(NA, c(dim(sim.forloop.result[[1]]$output),replications))
#   tfr <- tframe(sim.forloop.result[[1]]$output)
#   for (i in 1:replications) result[,,i] <- sim.forloop.result[[i]]$output
#  }
# else {
   #r <- simulate(model, list(...))$output
   r <- do.call("simulate", arglist)$output
   result <- array(NA, c(dim(r),replications))
   tfr <- tframe(r)
   result[,,1] <- r
   if (1 < replications)
     for (i in 2:replications) 
        result[,,i] <- outputData(do.call("simulate", arglist))
#  }
# default does not work on array seriesNames(result) <- seriesNamesOutput(model)
if (length(seriesNamesOutput(model)) != dim(result)[2])
 stop("length of names (",length(seriesNamesOutput(model)),
      ") does not match number of series(",dim(result)[2],").")
attr(result,"seriesNames") <- seriesNamesOutput(model)

result <- tframed(result, tfr)  # my more general multidimensional ts
invisible(classed( # constructor MonteCarloSimulations
         list(simulations=result, model=model, rng=rng, simulation.args=simulation.args,
              description = "data generated by MonteCarloSimulations.TSmodel"),
   c("MonteCarloSimulations") ))
}

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

print.MonteCarloSimulations <- function(x, digits=options()$digits, ...)
{#  (... further arguments, currently disregarded)
 cat("Simulation with RNG ", x$rng, " from model:\n")
 print(x$model)
 invisible(x)
}

nseriesOutput.MonteCarloSimulations <- function(x)
   {dim(x$simulations)[2]}

nseriesInput.MonteCarloSimulations <- function(x)
  {nseriesInput(x$simulation.args$data)}

tframe.MonteCarloSimulations <- function(x) tframe(x$simulations)

Tobs.MonteCarloSimulations <- function(x) Tobs(tframe(x))

seriesNamesOutput.MonteCarloSimulations <- function(x)
   {dimnames(x$simulations)[[2]]}

seriesNamesInput.MonteCarloSimulations <- function(x)
  {seriesNamesInput(x$simulation.args$data)}

testEqual.MonteCarloSimulations <- function(obj1, obj2, fuzz=1e-16)
 {if (length(obj1$result) != length(obj2$result)) r <- FALSE
  else  r <- all(fuzz > abs(obj1$simulations - obj2$simulations))
  r
 }


summary.MonteCarloSimulations <- function(object,
        series=NULL, periods=1:3, ...)
 {#  (... further arguments, currently disregarded)
  stats <- NULL
  if (!is.null(series))
    {if (dim(object$simulations)[3] <20) 
        warning("SD calculation is not very good with so few simulations.")
     names <- seriesNamesOutput(object)
     if(is.null(names)) names <- seriesNamesOutput(object$model)
     if (!is.numeric(series)) series <-match(series, names)
     names <- names[series]
     mn<-apply(object$simulations[periods,series,,drop=FALSE],c(1,2),mean)
     sd<-apply(object$simulations[periods,series,,drop=FALSE],c(1,2),var)^0.5
     stats <- rbind(mn,sd) 
     dimnames(stats)<- list(c(paste("mean period", periods), 
                          paste("S.D. period",periods)), names)
    }
  classed(list(  # constructor summary.MonteCarloSimulations
     description=object$description,
     sampleT=dim(object$simulations)[1],
     p=      dim(object$simulations)[2], 
     simulations=dim(object$simulations)[3], 
     summary.stats=stats,
     rng=getRNG(object)), 
  "summary.MonteCarloSimulations")
}


print.summary.MonteCarloSimulations <- function(x, digits=options()$digits, ...)
 {#  (... further arguments, currently disregarded)
  cat("Object class MonteCarloSimulations\n")
  cat(x$description, "\n")
  cat("Tobs=",x$sampleT, "variables=", x$p,"simulations=",x$simulations,"\n")
  cat("rng= ", x$rng, "\n")
  if (!is.null(x$summary.stats))   print(x$summary.stats, digits=digits)
  invisible(x)
 }


tfplot.MonteCarloSimulations <- function(x, 
    tf=tframe(x$simulations), start=tfstart(tf), end=tfend(tf),
    series=seq((dim(x$simulations)[2])), 
    select.simulations=seq(dim(x$simulations)[3]),
    graphs.per.page=5, mar=par()$mar, ...)
  {#  (... further arguments, currently disregarded)
   names <- seriesNames(x$simulations)
   tf.p <- tframe(x$simulations) # actual may be different (but not in default)
   Ngraph <- min(length(series), graphs.per.page)
   old.par <-par(mfcol = c(Ngraph, 1), mar= mar, no.readonly=TRUE) #c(5.1,6.1,4.1,2.1))
   on.exit(par(old.par))
   #zz<- matrix(NA, dim(sim)[1], length(x$simulations))
   if (!is.numeric(series)) series <- match(series, names)
   for(i in series) 
        {zz <- (x$simulations)[,i,select.simulations]
         tframe(zz) <- tf.p
	 seriesNames(zz) <- NULL #otherwise name length does not match in tfwindow.default(zz)
         tfOnePlot(zz,start=start,end=end, ylab=names[i]) #tsplot
         if(i == series[1])  title(main = "Monte Carlo Simulations")
        }
   invisible()
}



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

#distribution.TSdata <- function(obj, bandwidth=0.2, 
#        select.inputs = seq(length= nseriesInput(obj)),
#        select.outputs= seq(length=nseriesOutput(obj)), ...)
distribution.TSdata <- function(obj, ..., bandwidth=0.2, 
        select.inputs = seq(length= nseriesInput(obj)),
        select.outputs= seq(length=nseriesOutput(obj)))
  {#  (... further objects, currently disregarded)
   if (0 !=  nseriesInput(obj))
      distribution( inputData(obj), bandwidth=bandwidth, series=select.inputs)
   if (0 != nseriesOutput(obj))
      distribution(outputData(obj), bandwidth=bandwidth, series=select.outputs)
   invisible(obj)
  }

# this should be something like distribution.tframed if distribution becomes
#  generic in EstEval.  ALSO SEE distribution.factorsEstEval RE ...
#distribution.default <- function(obj, bandwidth=0.2, series=NULL, ...)
distribution.default <- function(obj, ..., bandwidth=0.2, series=NULL)
  {#  (... further objects, currently disregarded)
   # obj should be a ts matrix (perhaps this should be a tf method).
   # If series is NULL then all series are ploted.
   # note that this graphic can be fairly misleading:
   #    distribution(runif(1000))  should be uniform
   names <- seriesNames(obj)
   if (!is.matrix(obj) ) obj <- matrix(obj, length(obj), 1)
   if(!is.null(series)) 
     {obj   <-   obj[,series, drop=FALSE]
      names <- names[series]
     }
   par(mfcol=c(ncol(obj),1))
   for ( i in 1:ncol(obj))
      {if      (exists("density")) rd <- density(obj[,i], bw= bandwidth)
       #  'is.R' is deprecated in R 4.3.3 2024-02-02, replace with TRUE
       else if (exists("ksmooth") & !TRUE) rd <- ksmooth(obj[,i], bandwidth=bandwidth) 
       else     stop("Neither ksmooth nor density are available.")
       plot(rd, type="l", ylab="density", ylim=c(0, max(rd$y)), xlab=names[i] )
      }
   invisible()
  }


distribution.MonteCarloSimulations <- function(obj,
     series=seq(dim(obj$simulations)[2]),
     x.sections=TRUE, periods=1:3, graphs.per.page=5, ...)
  {#  (... further arguments, currently disregarded)
  if (dim(obj$simulations)[3] <20) 
     warning("This is not very good with so few simulations.")
   names <- seriesNamesOutput(obj)
   if(is.null(names)) names <- seriesNamesOutput(obj$model)
   if (!is.numeric(series)) series <- match(series, names)
   names <- names[series]
   Ngraph <- min(length(series), graphs.per.page)
   if (x.sections)
       {data <- obj$simulations[periods, series,,drop=FALSE]
        old.par <-par(mfrow =c(Ngraph, length(periods)),
	              mar=c(5.1,6.1,4.1,2.1), no.readonly=TRUE)
       }
   else 
       {old.par <-par(mfrow =c(Ngraph, 1), mar=c(5.1,6.1,4.1,2.1), no.readonly=TRUE)
        mn <- apply(obj$simulations[, series,,drop=FALSE], c(1,2), mean)
        sd <- apply(obj$simulations[, series,,drop=FALSE], c(1,2), var) ^ 0.5
        plt <- array(c(mn, mn+sd, mn-sd, mn+2*sd, mn-2*sd), c(dim(mn),5))
        tf.p <- tframe(obj$simulations)
      }
   on.exit(par(old.par))
   for (i in 1:length(series)) 
    {if (x.sections)
        {for (j in 1:length(periods)) 
          {if (exists("density")) plot(density(data[j,i,]), # -mean[j,i]
                 type="l",xlab=paste("Period",periods[j]),ylab=names[i])
           else if (exists("ksmooth") & !TRUE) plot(ksmooth(data[j,i,], # -mean[j,i],
                              bandwidth=var(data[j,i,])^0.5, kernel="parzen"),
                        type="l",xlab=paste("Period",periods[j]),ylab=names[i])
           else
        stop("Neither ksmooth nor density are available to calculate the plot.")
           if ((i == 1) & (j ==length(periods)%/%2))
              title(main = "kernel estimate of distributions")
           }
        }
     else
        {pl <-plt[,i,]
         tframe(pl) <- tf.p
         tfplot(pl, type="l", lty= c(1,3,3,2,2), ylab=names[i]) #tsplot
         if (i == 1) title(main = "Simulation mean, 1 & 2 S.D. estimates")
        }
    }
  invisible()
  }

############################################################################
#
#       methods for EstEval.  <<<<<<<<<<
#
############################################################################

#e.bb.ar.100 <- EstEval( mod2, replications=100, 
#               estimation.args=list(estimation="estVARXar", verbose=F))

#e.bb.ls.over <- EstEval( simple.mod, replications=100, 
#   estimation.args=list(estimation="estVARXls", max.lag=6, verbose=F), 
#   criterion="coef")


EstEval <- function( model, replications=100, rng=NULL, quiet=FALSE, 
                       simulation.args=NULL,
                       estimation=NULL, estimation.args=NULL, 
                       criterion ="coef", criterion.args =NULL) 
#		       Spawn=if (exists(".SPAWN")) .SPAWN else FALSE)
{
 if(is.null(estimation)) stop("estimation method must be specified.")
 if (is.EstEval(model) | is.MonteCarloSimulations(model))
   {rng  <- getRNG(model)
    model<- TSmodel(model)
   }
  
 truth <- do.call(criterion, append(list(model), criterion.args))

 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 (!Spawn) {
    if(!quiet) cat("Calculating ", replications, " estimates.\n")
    result <- vector("list",replications)
    for (i in 1:replications)
       {data <- do.call("simulate", append(list(model), simulation.args))
	m   <-  do.call(estimation, append(list(data),  estimation.args))
        result[[i]]<-do.call(criterion, append(list(m), criterion.args))
       }
#   }
# else {
#     if(!quiet)
#	 cat("Spawning processes to calculate ", replications, " estimates.\n")
#     est.forloop <- function(estimation, estimation.args, model, 
#			      simulation.args, criterion, criterion.args)
#	{data <- do.call("simulate", append(list(model), simulation.args))
#	 m   <-  do.call(estimation, append(list(data),  estimation.args))
#	 do.call(criterion, append(list(m), criterion.args))
#	}
#     assign("est.forloop", est.forloop, where = 1)
#     assign("est.forloop.n", replications, where = 1)
#     assign("est.forloop.result", list(NULL), where = 1)
#     assign("est.forloop.estimation", estimation, where = 1)
#     assign("est.forloop.model", model, where = 1)
#     assign("est.forloop.simulation.args", simulation.args, where = 1)
#     assign("est.forloop.criterion", criterion, where = 1)
#     assign("est.forloop.estimation.args", estimation.args, where = 1)
#     assign("est.forloop.criterion.args", criterion.args, where = 1)
#
#     on.exit(remove(c("est.forloop", "est.forloop.i", "est.forloop.n",
#	  "est.forloop.result",  "est.forloop.estimation","est.forloop.model",
#	  "est.forloop.simulation.args", "est.forloop.criterion", 
#	  "est.forloop.estimation.args","est.forloop.criterion.args"),where = 1))
#
#     For(est.forloop.i = 1:est.forloop.n, est.forloop.result[[est.forloop.i ]]<-
#	est.forloop(est.forloop.estimation, est.forloop.estimation.args, est.forloop.model, 
#	est.forloop.simulation.args, est.forloop.criterion, est.forloop.criterion.args),
#	first=options(warn=-1), sync = TRUE)
#     result<-est.forloop.result
#    }
invisible(classed( # constructor EstEval (EstEvals)
      list(result=result,truth=truth,model=model,
           rng=rng, version=version,
           estimation=estimation, estimation.args=estimation.args,
            criterion=criterion,   criterion.args=criterion.args, 
            simulation.args=simulation.args),
     c(paste(criterion,"EstEval",sep=""), "EstEval")))
}


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

testEqual.EstEval <- function(obj1, obj2, fuzz = 0)
 {all(as.character(obj1) == as.character(obj2))}

print.EstEval <- function(x, digits=options()$digits, ...)
{#  (... further arguments, currently disregarded)
 cat("Estimation evaluation with model:\n")
 print(x$model, digits=digits)
 cat("Evaluation criterion: ",x$criterion, "\n")
 invisible(x)
}

summary.EstEval <- function(object, ...)
 {#  (... further arguments, currently disregarded)
  classed(list( # constructor summary.EstEval
     class=class(object),
     estimation=object$estimation,
     estimation.args= if(!is.list((object$estimation.args)[[1]]))
                           object$estimation.args    else    NULL,
     criterion=object$criterion,
     labels=object$labels,
     criterion.args=object$criterion.args,
     replications=length(object$result), 
     true.model=object$model,
     rng=getRNG(object)), 
  "summary.EstEval")
 }


print.summary.EstEval <- function(x, digits=options()$digits, ...)
{ #  (... further arguments, currently disregarded)
  cat("Object of class: ", x$class, "\n")
  cat("Evaluation of `",x$estimation,"'")
  if(!is.list((x$estimation.args)[[1]]))
    {cat( " estimation with argument ") 
     cat(labels(x$estimation.args),"= `",x$estimation.args,"'")
    }
  cat("\n")
  cat("using criterion `", x$criterion, "' with argument ")
  cat(x$labels," = `", x$criterion.args, "'\n")
  cat(x$replications, " replications, RNG = ", x$rng, "\n")
  cat("true model:\n")
  print(x$model)
  invisible(x)
 }

tfplot.EstEval <- function(x, tf=NULL, start=tfstart(tf), end=tfend(tf),
        truth= if(is.TSdata(x$truth)) outputData(x$truth) else x$truth,
        series = seq(length=nseries(truth)),
	Title="Estimated (and true) results",
        ylab = seriesNames(truth), remove.mean=FALSE,
	graphs.per.page=5, mar=par()$mar, reset.screen=TRUE, ...)
 {#  (... further arguments, currently disregarded)
  N<-length(x$result)
  old.par <- par(mfcol = c(min(length(series), graphs.per.page), 1),
           no.readonly=TRUE)
  on.exit(par(old.par))
  tf <- tframe(truth)
  truth <- unclass(truth)
  if (remove.mean) 
    {truth <- sweep(truth,2, colMeans(truth))
     ylab <- paste(ylab, "(mean removed)")
    }
  r <- matrix(NA,Tobs(truth), N)
  for (j in series) {
  	for (i in 1:N)  r[,i] <- x$result[[i]][,j]
	if (remove.mean) r <- sweep(r,2, colMeans(r))
  	tfOnePlot(tframed(tbind(truth[,j], r), tf), ylab= ylab[j],
	              lty=c(1,rep(2, N)), col=c("black",rep("red",N)),
		      start=start,end=end)
        if(!is.null(Title) && (j == series[1]) && (is.null(options()$PlotTitles)
                || options()$PlotTitles)) title(main = Title)
        }
  invisible()
 }



############################################################################
#
#       methods for rootsEstEval  (EstEval)  <<<<<<<<<<
#
############################################################################

summary.rootsEstEval <- function(object, verbose=TRUE, ...)
{#  (... further arguments, currently disregarded)
  nxt <- if (verbose) NextMethod("summary") else NULL
  if (! verbose) conv <- NULL
  else 
    {if (!is.null(object$result.conv)) conv <- NULL
     else                              conv <- object$result.conv
    }
  N <- length(object$result)
  p <- 0
  for (i in 1:N) p <- max(p, length((object$result)[[i]]))
  r <- matrix(NA, N, p)
  for (i in 1:N) r[i,1:length((object$result)[[i]])] <- (object$result)[[i]]
  m <- colSums(r)/N
  cov <- r- t(matrix(object$truth, p, N))
  cov <- (t(Conj(cov)) %*% cov)/(N-1)
  ecov <- r- t(matrix(m, p, N))
  ecov <- (t(Conj(ecov)) %*% ecov)/(N-1)
  classed(list( # constructor summary.summary.rootsEstEval
     nxt=nxt,
     conv=conv,
     true.criterion=object$truth,
     mean=m,
     cov=cov,
     ecov=ecov),
  "summary.summary.rootsEstEval")
 }

print.summary.rootsEstEval <- function(x, digits=options()$digits, ...)
 {#  (... further arguments, currently disregarded)
  if (!is.null(x$nxt)) print(x$nxt)
  if (!is.null(x$conv))
        {if(all(x$conv)) cat("All estimates converged.\n")
         else cat(sum(!x$conv)," estimates did not converge!\n")
        }
  cat("\nTrue model criterion mean: ",x$true.criterion,"\n")
  cat("Sampling estimate of mean: ",x$mean,"\n")
  cat("Estimate of sampling covariance [e*Conj(t(e))] using true model:\n")
  print(x$cov)
  cat("\nEstimate of sampling covariance (using sample mean and not the true model):\n")
  print(x$ecov)
  invisible(x)
 }


tfplot.rootsEstEval <- function(x, ...)UseMethod("plot.rootsEstEval")

plot.rootsEstEval <- function(x, complex.plane=TRUE, cumulate=TRUE, norm=FALSE,
     bounds=TRUE, transform=NULL, invert=FALSE, Sort=TRUE, ...)
{#  (... further arguments, currently disregarded)
   N<-length(x$result)
   n <- 0
   for (i in 1:N) n <- max(n, length((x$result)[[i]]))
   r <- matrix(0,N, n) 
   for (i in 1:N) r[i,1:length((x$result)[[i]])] <- (x$result)[[i]]
   true.lines <- c(x$truth, rep(0,n-length(x$truth)))
   if (invert)
         {true.lines <- 1/true.lines
          r <- 1/r
         }
   if(!is.null(transform)) 
         {r <- do.call(transform,list(r))
          true.lines <-do.call(transform,list(true.lines))
         }
   if (complex.plane)
    {#plot.roots(x$truth, pch="o")
     plot(x$truth, pch="o")
     for (i in 1:N) addPlotRoots(r[i,], pch="*") 
     addPlotRoots(0, pch="+") # addPlotRoots(0+0i, pch="+")
    }
  else
     {if (Sort)
        {r <- t(apply(r,1,sort))
         true.lines <- sort(true.lines)
        }
      if (cumulate) r <- apply(r,2,cumsum)/matrix(1:N,N,ncol(r))
      else     r <- r 
      if(is.complex(r))
         { r <- cbind(Re(r), Im(r))
          true.lines <-c(Re(true.lines),Im(true.lines))
         }
      r[is.infinite(r)] <- 0
      true.lines <-t(matrix(true.lines, length(true.lines),N))
      matplot(x=seq(nrow(r)), y=cbind(0,true.lines, r), type="l",
             lty=c(1,rep(3,dim(true.lines)[2]), rep(2,dim(r)[2])) )
     }
  invisible(r)
}

roots.rootsEstEval <- function(obj, ...)   {obj}

distribution.rootsEstEval <- function(obj, ..., mod=TRUE, invert=FALSE, Sort=FALSE, 
    bandwidth=0.2, select=NULL)
{#  (... further objectss, currently disregarded)
 # if mod is true the modulus is used, otherwise real and imaginary are separated.
 # if invert is true the reciprical is used.
 # if Sort is true then sort is applied (before cumulate). This is of particular interest
 #   with estimation methods like black.box which may not return parameters
 #   of the same length or in the same order.
 # If select is not NULL then only the indicated roots are plotted. 
 #     ie - select=c(1,2)  will plot only the two largest roots
      N<-length(obj$result)
      n <- 0
      for (i in 1:N) n <- max(n, length((obj$result)[[i]]))
      r <- matrix(0,N,n)
      for (i in 1:N) r[i,] <- c((obj$result)[[i]], 
                                rep(0,n-length((obj$result)[[i]])))
      true.lines <- c(obj$truth, rep(0,n-length(obj$truth)))
      if (invert)
         {true.lines <- 1/true.lines
          r <- 1/r
         }
      if(mod) 
         {r <- Mod(r) 
          true.lines <-Mod(true.lines)
          xlab <-"Mod root "
         }
      else
         { r <- cbind(Re(r), Im(r))
          true.lines <-c(Re(true.lines),Im(true.lines))
          xlab <-"Real part root "
         }
      r[is.infinite(r)] <- 0
      if (Sort)
        {r <- t(apply(r,1,sort))
         true.lines <- sort(true.lines)
        }
      if(!is.null(select)) r <- r[,select, drop=FALSE]
      par(mfcol=c(dim(r)[2],1))
      for ( i in 1:dim(r)[2])
         {if      (exists("ksmooth") & !TRUE) rd <- ksmooth(r[,i], bandwidth=bandwidth) 
          else if (exists("density")) rd <- density(r[,i], bw= bandwidth)
          else
        stop("Neither ksmooth nor density are available to calculate the plot.")
          if (i > n) xlab <-"Imaginary part root "
          plot(rd, type="l", ylab="density", ylim=c(0, max(rd$y)),
               xlab=paste(xlab, n-(-i%%n)) )
          lines(rep(true.lines[i],2),c(1,0))
         }
      invisible()
}


############################################################################
#
#       methods for coefEstEval (EstEval)  <<<<<<<<<<
#
############################################################################

summary.coefEstEval <- function(object, verbose=TRUE, ...)
  {classed(summary.rootsEstEval(object, verbose=verbose), "summary.coefEstEval")} # constructor
#  (... further arguments, currently disregarded)


print.summary.coefEstEval <- function(x, digits=options()$digits, ...)
  UseMethod("print.summary.rootsEstEval")
#  (... further arguments, currently disregarded)


tfplot.coefEstEval <- function(x, cumulate=TRUE, norm=FALSE, bounds=TRUE,
        invert=FALSE, Sort=FALSE, graphs.per.page = 5, ...){
      #  (... further arguments, currently disregarded)
      N<-length(x$result)
      n <- 0
      for (i in 1:N) n <- max(n, length((x$result)[[i]]))
      r <- matrix(0,N,n)
      plottrue <- TRUE
      for (i in 1:N) {
	ni <- length((x$result)[[i]])
	r[i,1:ni] <- (x$result)[[i]]
	if (ni != n) plottrue <- FALSE
	}
      if (invert) r <- 1/r
      if(norm)    r <- matrix((rowSums(r^2))^.5, N,1)
      r[is.infinite(r)] <- 0
      if (Sort) r <- t(apply(r,1,sort))
      if (n != length(x$truth)) plottrue <- FALSE
      if (plottrue) {
	true.lines <- c(x$truth)
	if (invert) true.lines <- 1/true.lines
	if(norm) true.lines <-sum(true.lines^2)^.5
	if (Sort) true.lines <- sort(true.lines)
	true.lines <-t(matrix(true.lines, length(true.lines),N))
	if (bounds) {
		z  <- r-true.lines
		#Om <- t(z) %*% z/(nrow(z)-1)
		Om <- crossprod(z, z)/(nrow(z)-1)
		Om <- diag(Om)^.5
		Om <- t(matrix(Om, length(Om), N))
		Om <- Om/matrix((1:N)^.5 , N, ncol(Om))
        }
      }
      if (cumulate) r<- apply(r,2,cumsum)/matrix(1:N,N,ncol(r))
      seriesNames(r) <- paste("parm", seq(ncol(r)))
#      matplot(x=matrix(seq(nrow(r)),nrow(r),1), y=cbind(0,true.lines,r, Om), 
#              type="l", lty=c(1,rep(3,dim(true.lines)[2]), rep(4,dim(r)[2]), 
#                     rep(2,2*dim(r)[2])))
      if (plottrue & bounds)
         tfplot(r, true.lines,true.lines + Om, true.lines - Om,
            graphs.per.page = graphs.per.page)
      else if (plottrue) tfplot(r, true.lines, graphs.per.page=graphs.per.page)
      else               tfplot(r, graphs.per.page = graphs.per.page)
      invisible(r)
}

roots.coefEstEval <- function(obj, criterion.args=NULL, ...)
{# extract roots criterion 
  model <- obj$model
  truth <-do.call("roots", append(list(model), criterion.args))
  r <- NULL
  for (m in obj$result)
    {coef(model) <- m
     model <- setArrays(model)
     r <- append(r, 
           list(do.call("roots", append(list(model), criterion.args))))
    }
  ok <- TRUE
  for (m in obj$result)
     ok <- ok & (length(coef(model)) == length(m))  # not perfect but ...
  if (!ok) warning("Parameters do not all correspond to given true model.")
  obj$result<-r
  obj$truth <-truth
  obj$criterion<-"roots"
  obj$criterion.args <-criterion.args
  invisible(classed(obj, c("rootsEstEval","EstEval")))# constructor
}

distribution.coefEstEval <- function(obj, ...,  Sort=FALSE, bandwidth=0.2,
	graphs.per.page=5)
{#  (... further objects, currently disregarded)
 # if Sort is true then sort is applied (before ave). This is of particular interest
 #   with estimation methods like black.box which may not return parameters
 #   of the same length or in the same order.
      N<-length(obj$result)
      n <- length(obj$truth)
      for (i in 1:N) n <- max(n, length((obj$result)[[i]]))
      if (n == length(obj$truth)) plottrue <- TRUE
      else {
         warning("Number of true parameters does not match number estimated.")
	 plottrue <- FALSE
      }
      r <- matrix(0,N,n)
      for (i in 1:N) r[i,1:length((obj$result)[[i]])] <- (obj$result)[[i]]
      true.lines <- c(obj$truth, rep(0,n-length(obj$truth)))
      if (Sort)
        {r <- t(apply(r,1,sort))
         true.lines <- sort(true.lines)
        }
      xlab <-"parameter "
      old.par <- par(mfcol=c(min(graphs.per.page, ncol(r)),1),
                    mar=c(5.1, 4.1, 4.1, 2.1), no.readonly=TRUE)
      on.exit(par(old.par))
      for ( i in 1:ncol(r))
         {if (TRUE)     rd <- density(r[,i], bw=bandwidth)
          else            rd <- ksmooth(r[,i], bandwidth=bandwidth) # Splus
          plot(rd, type="l", ylim=c(0, max(rd$y)),
               ylab="density",  xlab=paste(xlab, i) , main="")
          if (plottrue) lines(rep(true.lines[i],2),c(1,0))
         }
      invisible()
}

TSmodel.coefEstEval <- function(obj, ...)
{# rebuild model from coef
  model <- obj$model
  truth <-TSmodel(model)
  r <- NULL
  for (m in obj$result)
    {coef(model) <- m
     model <- setArrays(model)
     r <- append(r, list(model))
    }
  ok <- TRUE
  for (m in obj$result)
     ok <- ok & (length(coef(model)) == length(m))  # not perfect but ...
  if (!ok) warning("Parameters do not all correspond to given true model.")
  obj$result<-r
  obj$truth <-truth
  obj$criterion<-"TSmodel"
  obj$criterion.args <-NULL
  invisible(classed(obj, c("TSmodelEstEval","EstEval")))
}

TSestModel.coefEstEval <- function(obj)
{# rebuild ... 
  model <- obj$model
  truth <-l(TSmodel(model), data)   # need to regenerate data
  r <- NULL
  for (m in obj$result)
    {coef(model) <- m
     model <- l( setArrays(model), data)
     r <- append(r, list(model))
    }
  ok <- TRUE
  for (m in obj$result)
     ok <- ok & (length(coef(model)) == length(m))  # not perfect but ...
  if (!ok) warning("Parameters do not all correspond to given true model.")
  obj$result<-r
  obj$truth <-truth
  obj$criterion<-"TSestModel"
  obj$criterion.args <-NULL
  invisible(classed(obj, c("TSestModelEstEval","EstEval")))
}

############################################################################
#
#       methods for TSmodelEstEval (EstEval)  <<<<<<<<<<
#
############################################################################

summary.TSmodelEstEval <- function(object, ...)
 {#  (... further arguments, currently disregarded)
  if (is.null((object$result)[[1]]$converged)) conv <- NULL
  else
    {conv <- rep(NA,length(object$result))
     for (i in 1:length(conv)) conv[i] <- (object$result)[[i]]$converged
    }
 # summary(coef(object))
 #summary(roots(object))  these are slow

  classed(list( # constructor summary.TSmodelEstEval
     class=class(object),
     conv=conv,
     default=summary.default(object)),
  "summary.TSmodelEstEval")
 }

print.summary.TSmodelEstEval <- function(x, digits=options()$digits, ...)
{#  (... further arguments, currently disregarded)
  cat("Object of class: ",x$class, "\n")
  if (!is.null(x$conv))
        {if(all(x$conv)) cat("All estimates converged.\n")
         else cat(sum(!x$conv)," estimates did not converge!\n")
        }
  print(x$default)
  invisible(x)
}


coef.TSmodelEstEval <- function(object, criterion.args=NULL, ...)
{#  (... further arguments, currently disregarded)
 # extract parameters from models in the list and 
 #   return a list of class coefEstEval EstEval 
 # criterion.args is not used. It is provided only so calls from 
 #   summary.TSmodelEstEval can provide this argument.
  truth <-coef(object$truth)
  r <- NULL
  for (m in object$result) 
     r <- append(r,list(coef(m)))
  if (! is.null((object$result)[[1]]$converged))
    {rc <- rep(NA,length(object$result))
     for (i in 1:length(rc)) rc[i] <- (object$result)[[i]]$converged
     object$result.conv<-rc
    }
  object$result<-r
  object$truth <-truth
  object$criterion<-"coef"
  object$criterion.args <-criterion.args
  invisible(classed(object, c("coefEstEval","EstEval")))
}

roots.TSmodelEstEval <- function(obj, criterion.args=list( randomize=TRUE), ...)
{# extract roots criterion 
  truth <-do.call("roots", append(list(obj$truth), criterion.args))
  r <- NULL
  for (m in obj$result)
     r <- append(r, 
           list(do.call("roots", append(list(m), criterion.args))))
  if (! is.null((obj$result)[[1]]$converged))
    {rc <- rep(NA,length(obj$result))
     for (i in 1:length(rc)) rc[i] <- (obj$result)[[i]]$converged
     obj$result.conv<-rc
    }
  obj$result<-r
  obj$truth <-truth
  obj$criterion<-"roots"
  obj$criterion.args <-criterion.args
  invisible(classed(obj, c("rootsEstEval","EstEval")))
}


tfplot.TSmodelEstEval <- function(x, graph.args=NULL,
                       criterion ="coef", criterion.args=NULL, ...){
  #  (... further arguments, currently disregarded)
  # extract criterion and pass to another method with graph.args
  r <- do.call(paste(criterion,".TSmodelEstEval", sep=""), 
               append(list(x), list(criterion.args=criterion.args)))
  do.call("tfplot", append(list(r), graph.args))
  invisible(r)
  }

############################################################################
#
#       methods for TSestModelEstEval (EstEval)   <<<<<<<<<<
#
############################################################################

summary.TSestModelEstEval <- function(object, ...)
  {classed(summary.TSmodelEstEval(object), "summary.TSestModelEstEval") }  # constructor 
#  (... further arguments, currently disregarded)

print.summary.TSestModelEstEval <- function(x, digits=options()$digits, ...)
   UseMethod("print.summary.TSmodelEstEval")
#  (... further arguments, currently disregarded)


coef.TSestModelEstEval <- function(object, criterion.args=NULL, ...)
{#  (... further arguments, currently disregarded)
 # extract parameters from models in the list and convergence info.
 #   return a list of class coefEstEval EstEval
 # criterion.args is not used. It is provided only so calls from 
 #   summary.TSmodelEstEval can provide this argument.
  truth <-coef(object$truth)
  r <- NULL
  for (m in object$result) r <- append(r,list(coef(m)))
  rc <- rep(NA,length(object$result))
  for (i in 1:length(rc)) rc[i] <- (object$result)[[i]]$converged
  object$result<-r
  object$result.conv<-rc
  object$truth <-truth
  object$criterion<-"coef"
  object$criterion.args <-criterion.args
  invisible(classed(object, c("coefEstEval","EstEval")))
}

roots.TSestModelEstEval <- function(obj, criterion.args=NULL, ...)
{# extract roots criterion 
  truth <-do.call("roots", append(list(obj$truth), criterion.args))
  r <- NULL
  for (m in obj$result)
     r <- append(r, 
             list(do.call("roots", append(list(m), criterion.args))))
  rc <- rep(NA,length(obj$result))
  for (i in 1:length(rc)) rc[i] <- (obj$result)[[i]]$converged
  obj$result<-r
  obj$result.conv<-rc
  obj$truth <-truth
  obj$criterion<-"roots"
  obj$criterion.args <-criterion.args
  invisible(classed(obj, c("rootsEstEval","EstEval")))
}


tfplot.TSestModelEstEval <- function(x, graph.args=NULL,
                       criterion ="coef", criterion.args=NULL, ...){
  #  (... further arguments, currently disregarded)
  # extract criterion and pass to another method with graph.args
  r <- do.call(paste(criterion,".TSestModelEstEval", sep=""), 
               append(list(x), list(criterion.args=criterion.args)))
  do.call("tfplot", append(list(r), graph.args))
  invisible(r)
  }


############################################################################
#
#  methods for generating test data
#
############################################################################

genMineData <- function(umodel, ymodel, uinput=NULL, sampleT=100, 
	unoise=NULL, usd=1,ynoise=NULL, ysd=1, rng=NULL)
{if (is.TSestModel(umodel)) umodel <- TSmodel(umodel)
 if (is.TSestModel(ymodel)) ymodel <- TSmodel(ymodel)
 if(!is.TSmodel(umodel)) stop("genMineData expecting a TSmodel.")
 if(!is.TSmodel(ymodel)) stop("genMineData expecting a TSmodel.")

 if (nseriesInput(ymodel) != nseriesOutput(umodel))
   stop("umodel output dimension must equal ymodel input dimension.")
 
 if(is.null(rng)) rng <- setRNG() # returns setting so don't skip if NULL
 else        {old.rng <- setRNG(rng);  on.exit(setRNG(old.rng))  }
 
 input <- outputData(simulate(umodel, input=uinput, sampleT=sampleT, 
                   noise=unoise, sd=usd))
 r <- TSdata(input  = input,
             output = outputData(simulate(ymodel, input=input,
	                             sampleT=sampleT, noise=ynoise, sd=ysd)) )
 r$umodel  <- umodel
 r$ymodel  <- ymodel
 r$uinput  <- uinput
 r$sampleT <- sampleT 
 r$unoise  <- unoise
 r$usd     <- usd
 r$ynoise  <- ynoise
 r$ysd     <- ysd
 r$rng     <- rng
 r
}

build.input.models <- function(data, max.lag=NULL)
{# make a list of univariate models, one for each series in inputData(data)
 #   for use by build.diagonal.model. 
 n <- nseriesInput(data)
 multi.models <- vector("list", n)
 for (i in seq(n))
   {d <-trimNA(TSdata(output= inputData(data, series=i)))
    multi.models[[i]] <- TSmodel(estVARXls(d, max.lag=max.lag))
   }
 multi.models
}

build.diagonal.model <- function(multi.models)
{# build one diagonal model from a list of models as returned  by 
 # build.input.models. Uses the AR part only. This can be used by genMineData.
 n <- length(multi.models)
 lag <- 0
 for (i in seq(n)) lag <- max(lag, dim(multi.models[[i]]$A)[1])
 p <- 0
 for (i in seq(n))  p  <- p + dim(multi.models[[i]]$A)[3]
 model <- array(0, c(lag, p,p))
 p <- 0
 for (i in seq(n))
   {pi <- dim(multi.models[[i]]$A)[3]
    li <- dim(multi.models[[i]]$A)[1]
    model[ 1:li, (p+1):(p+pi), (p+1):(p+pi)] <-  multi.models[[i]]$A
    p <- p + pi
   }
 ARMA(A= model, B=array(diag(1, p), c(1,p,p)))
}