Nothing
############################################################################
# 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)
else if (exists("ksmooth") & !is.R()) 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") & !is.R()) 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") & !is.R()) 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 (is.R()) 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)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.