ts.R

#  File src/library/stats/R/ts.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1995-2023 The R Core Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

start	  <- function(x, ...) UseMethod("start")
end	  <- function(x, ...) UseMethod("end")
frequency <- function(x, ...) UseMethod("frequency")
time	  <- function(x, ...) UseMethod("time")
window	  <- function(x, ...) UseMethod("window")
cycle     <- function(x, ...) UseMethod("cycle")
deltat    <- function(x, ...) UseMethod("deltat")

ts <- function(data = NA, start = 1, end = numeric(), frequency = 1,
	       deltat = 1, ts.eps  =  getOption("ts.eps"),
	       class = if(nseries > 1) c("mts", "ts", "matrix", "array") else "ts",
               names = if(!is.null(dimnames(data))) colnames(data)
                       else paste("Series", seq(nseries))
               )
{
    if(is.data.frame(data)) data <- data.matrix(data)
#   if(!is.numeric(data)) stop("'data'  must be a numeric vector or matrix")
    if(is.matrix(data)) {
	nseries <- ncol(data)
	ndata <- nrow(data)
        dimnames(data) <- list(NULL, names)
    } else {
	nseries <- 1
	ndata <- length(data)
    }
    if(ndata == 0) stop("'ts' object must have one or more observations")

    if(missing(frequency)) frequency <- 1/deltat
    else if(missing(deltat)) deltat <- 1/frequency

    if(frequency > 1 && 0 < (d <- abs(frequency - round(frequency))) && d < ts.eps)
	frequency <- round(frequency)

    ## FIXME: as.numeric
    if(!missing(start))
        start <- as.numeric(start)
    if(length(start) > 1L) {
## strange: this never checked for < 1!  commented for 1.7.0
##	if(start[2L] > frequency) stop("invalid start")
	start <- start[1L] + (start[2L] - 1)/frequency
    }
    ## FIXME: as.numeric
    if(!missing(end))
        end <- as.numeric(end)
    if(length(end) > 1L) {
##	if(end[2L] > frequency) stop("invalid end")
	end <- end[1L] + (end[2L] - 1)/frequency
    }
    if(missing(end))
	end <- start + (ndata - 1)/frequency
    else if(missing(start))
	start <- end - (ndata - 1)/frequency

    if(start > end) stop("'start' cannot be after 'end'")

    ## FIXME: as.numeric
    cycles <- as.numeric((end - start)*frequency) # as.n*(): get rid of date/time
    if(abs(round(cycles) - cycles) > ts.eps * max(cycles, 1))
    	stop("'end' must be a whole number of cycles after 'start'")
    nobs <- floor(cycles + 1.01)

    if(nobs != ndata)
	data <-
	    if(NCOL(data) == 1) {
		if(ndata < nobs) rep_len(data, nobs)
		else if(ndata > nobs) data[1L:nobs]
	    } else {
		if(ndata < nobs) data[rep_len(1L:ndata, nobs), ]
		else if(ndata > nobs) data[1L:nobs, ]
	    }
    ## FIXME: The following "attr<-"() calls C tspgets() which uses a
    ##  	fixed equivalent of ts.eps := 1e-5
    attr(data, "tsp") <- c(start, end, frequency) #-- order is fixed
    if(!is.null(class) && class[[1]] != "none") attr(data, "class") <- class
    ## if you alter the return structure, you also need to alter
    ## newBasic in methods/R/RClassUtils.R.  So please don't.
    data
}

tsp <- function(x) attr(x, "tsp")

`tsp<-` <- function(x, value)
{
    cl <- oldClass(x)
    attr(x, "tsp") <- value # does error-checking internally
    if (is.null(value)) {
        if (inherits(x, "ts"))
	    cl <- cl["ts" != cl]
        if (inherits(x, "mts"))
	    cl <- cl["mts" != cl]
        class(x) <- cl
    }
    x
}

hasTsp <- function(x)
{
    if(is.null(attr(x, "tsp")))
        attr(x, "tsp") <- c(1, NROW(x), 1)
    x
}

is.ts <- function(x) inherits(x, "ts") && length(x)

as.ts <- function(x, ...) UseMethod("as.ts")

as.ts.default <- function(x, ...)
{
    if (is.ts(x)) x
    else if(!is.null(xtsp <- tsp(x))) ts(x, xtsp[1L], xtsp[2L], xtsp[3L])
    else ts(x)
}

.cbind.ts <- function(sers, nmsers, dframe = FALSE, union = TRUE,
                      ts.eps = getOption("ts.eps"))
{
    nulls <- vapply(sers, is.null, NA)
    sers <- sers[!nulls]
    nser <- length(sers)
    if(nser == 0L) return(NULL)
    if(nser == 1L)
        return(if(dframe) as.data.frame(sers[[1L]]) else sers[[1L]])
    tsser <- vapply(sers, function(x) length(tsp(x)) > 0L, NA)
    if(!any(tsser))
        stop("no time series supplied")
    sers <- lapply(sers, as.ts)
    tsps <- sapply(sers[tsser], tsp)
    freq <- mean(tsps[3,])
    if(max(abs(tsps[3,] - freq)) > ts.eps)
        stop("not all series have the same frequency")

    ## cos(2pi ph) + 1i*sin(2pi ph); ph := phases
    eph <- exp(2i*(pi * apply(tsps, 2L, function(tsp) (tsp[1L]*tsp[3L]) %% 1)))
    if(max(Mod(eph - mean(eph))) > ts.eps)
    	stop("not all series have the same phase")
    if(union) {
        st <- min(tsps[1,])
        en <- max(tsps[2,])
    } else {
        st <- max(tsps[1,])
        en <- min(tsps[2,])
        if(st > en) {
            warning("non-intersecting series")
            return(NULL)
        }
    }
    p <- c(st, en, freq)
    n <- round(freq * (en - st) + 1)
    if(any(!tsser)) {
        ln <- vapply(sers[!tsser], NROW, 1)
        if(any(ln != 1L & ln != n)) # vectors, so not &&
            stop("non-time series not of the correct length")
        for(i in (1L:nser)[!tsser]) {
            sers[[i]] <- ts(sers[[i]], start=st, end=en, frequency=freq)
        }
        tsps <- sapply(sers, tsp)
    }
    nsers <- vapply(sers, NCOL, 1)
    if(dframe) {
	x <- setNames(vector("list", nser), nmsers)
    } else {
        ns <- sum(nsers)
        x <- matrix(, n, ns)
        cs <- c(0, cumsum(nsers))
        nm <- character(ns)
        for(i in 1L:nser)
            if(nsers[i] > 1) {
                cn <- colnames(sers[[i]]) %||% 1L:nsers[i]
                nm[(1+cs[i]):cs[i+1]] <- paste(nmsers[i], cn, sep=".")
            } else nm[cs[i+1]] <- nmsers[i]
        dimnames(x) <- list(NULL, nm)
    }
    for(i in 1L:nser) {
        if(union) {
            xx <-
                if(nsers[i] > 1)
                    rbind(matrix(NA, round(freq * (tsps[1,i] - st)), nsers[i]),
                          sers[[i]],
                          matrix(NA, round(freq * (en - tsps[2,i])), nsers[i]))
                else
                    c(rep.int(NA, round(freq * (tsps[1,i] - st))), sers[[i]],
                      rep.int(NA, round(freq * (en - tsps[2,i]))))
        } else {
            xx <- window(sers[[i]], st, en)
        }
        if(dframe) x[[i]] <- structure(xx, tsp=p, class="ts")
        else x[, (1+cs[i]):cs[i+1]] <- xx
    }
    if(dframe) as.data.frame(x)
    else ts(x, start=st, frequency=freq)
}

.makeNamesTs <- function(...)
{
    l <- as.list(substitute(list(...)))[-1L]
    nm <- names(l)
    fixup <- if(is.null(nm)) seq_along(l) else nm == ""
    ## <NOTE>
    dep <- sapply(l[fixup], function(x) deparse(x)[1L])
    ## We could add support for 'deparse.level' here by creating dep
    ## as in list.names() inside table().  But there is a catch: we
    ## need deparse.level = 2 to get the 'usual' deparsing when the
    ## method is invoked by the generic ...
    ## </NOTE>
    if(is.null(nm)) return(dep)
    if(any(fixup)) nm[fixup] <- dep
    nm
}

Ops.ts <- function(e1, e2)
{
    if(missing(e2)) {
        ## univariate operator
        NextMethod(.Generic)
    } else if(any(!nzchar(.Method))) {
        ## one operand is not a ts
        NextMethod(.Generic)
    } else {
        nc1 <- NCOL(e1)
        nc2 <- NCOL(e2)
        ## align e1 and e2
        e12 <- .cbind.ts(list(e1, e2),
                         c(deparse(substitute(e1))[1L],
                           deparse(substitute(e2))[1L]),
                         union = FALSE)
        e1 <- if(is.matrix(e1)) e12[,        1L:nc1 , drop = FALSE] else e12[, 1]
        e2 <- if(is.matrix(e2)) e12[, nc1 + (1L:nc2), drop = FALSE] else e12[, nc1 + 1]
        NextMethod(.Generic)
    }
}

cbind.ts <- function(..., deparse.level = 1) {
    if(deparse.level != 1) .NotYetUsed("deparse.level != 1")
    .cbind.ts(list(...), .makeNamesTs(...), dframe = FALSE, union = TRUE)
}

ts.union <- function(..., dframe = FALSE)
    .cbind.ts(list(...), .makeNamesTs(...), dframe = dframe, union = TRUE)

ts.intersect <- function(..., dframe = FALSE)
    .cbind.ts(list(...), .makeNamesTs(...), dframe = dframe, union = FALSE)

diff.ts <- function (x, lag = 1, differences = 1, ...)
{
    if (lag < 1 || differences < 1)
        stop("bad value for 'lag' or 'differences'")
    if (lag * differences >= NROW(x)) return(x[0L])
    ## <FIXME>
    ## lag() and its default method are defined in package ts, so we
    ## need to provide our own implementation.
    tsLag <- function(x, k = 1) {
        p <- tsp(x)
        tsp(x) <- p - (k/p[3L]) * c(1, 1, 0)
        x
    }
    r <- x
    for (i in 1L:differences) {
        r <- r - tsLag(r, -lag)
    }
    xtsp <- attr(x, "tsp")
    if(is.matrix(x)) colnames(r) <- colnames(x)
    ts(r, end = xtsp[2L], frequency = xtsp[3L])
}

na.omit.ts <- function(object, ...)
{
    tm <- time(object)
    xfreq <- frequency(object)
    ## drop initial and final NAs
    if(is.matrix(object))
        good <- which(apply(!is.na(object), 1L, all))
    else  good <- which(!is.na(object))
    if(!length(good)) stop("all times contain an NA")
    omit <- integer()
    n <- NROW(object)
    st <- min(good)
    if(st > 1) omit <- c(omit, 1L:(st-1))
    en <- max(good)
    if(en < n) omit <- c(omit, (en+1):n)
    cl <- attr(object, "class")
    if(length(omit)) {
        object <- if(is.matrix(object)) object[st:en,] else object[st:en]
        attr(omit, "class") <- "omit"
        attr(object, "na.action") <- omit
        tsp(object) <- c(tm[st], tm[en], xfreq)
        if(!is.null(cl)) class(object) <- cl
    }
    if(anyNA(object)) stop("time series contains internal NAs")
    object
}

is.mts <- function (x) is.ts(x) && inherits(x, "mts") && is.matrix(x)

start.default <- function(x, ts.eps = getOption("ts.eps"), ...)
{
    tsp <- attr(hasTsp(x), "tsp")
    is <- tsp[1L]*tsp[3L]
    if(abs(tsp[3L] - round(tsp[3L])) < ts.eps &&
       abs(is - round(is)) < ts.eps) {
	is <- floor(tsp[1L]+ts.eps)
	fs <- floor(tsp[3L]*(tsp[1L] - is)+0.001)
	c(is,fs+1)
    }
    else tsp[1L]
}

end.default <- function(x, ts.eps = getOption("ts.eps"), ...)
{
    tsp <- attr(hasTsp(x), "tsp")
    is <- tsp[2L]*tsp[3L]
    if(abs(tsp[3L] - round(tsp[3L])) < ts.eps &&
       abs(is - round(is)) < ts.eps) {
	is <- floor(tsp[2L]+ts.eps)
	fs <- floor(tsp[3L]*(tsp[2L] - is)+0.001)
	c(is, fs+1)
    }
    else tsp[2L]
}

frequency.default <- function(x, ...)
    if(!is.null(xtsp <- attr(x, "tsp"))) xtsp[3L] else 1

deltat.default <- function(x, ...)
    if(!is.null(xtsp <- attr(x, "tsp"))) 1/xtsp[3L] else 1

time.default <- function (x, offset = 0, ts.eps = getOption("ts.eps"), ...)
{
    xtsp <- attr(hasTsp(x), "tsp")
    y <- seq.int(xtsp[1L], xtsp[2L], length.out = NROW(x)) + offset/xtsp[3L]
    if(ts.eps > 0) {
        iy <- round(y)
        nearI <- abs(iy - y) < ts.eps
        y[nearI] <- iy[nearI]
    }
    tsp(y) <- xtsp
    y
}

time.ts <- function (x, ...) as.ts(time.default(x, ...))

cycle.default <- function(x, ...)
{
    p <- tsp(hasTsp(x))
    m <- round((p[1L] %% 1) * p[3L])
    x <- (1L:NROW(x) + m - 1) %% p[3L] + 1
    tsp(x) <- p
    x
}

cycle.ts <- function (x, ...) as.ts(cycle.default(x, ...))

print.ts <- function(x, calendar, ...)
{
    x <- as.ts(x)
    ## sanity check
    Tsp <- tsp(x)
    if(is.null(Tsp)) {
	warning("series is corrupt, with no 'tsp' attribute")
	print(unclass(x), ...)
	return(invisible(x))
    }
    nn <- 1 + round((Tsp[2L] - Tsp[1L]) * Tsp[3L])
    if(NROW(x) != nn) {
        warning(gettextf("series is corrupt: length %d with 'tsp' implying %d",
                         NROW(x), nn), domain=NA, call.=FALSE)
        calendar <- FALSE
    }
    fr.x <- frequency(x)
    if(missing(calendar))
	calendar <- any(fr.x == c(4,12)) && length(start(x)) == 2L
    if(!calendar) {
        if(fr.x != 1)
            cat("Time Series:\nStart =", deparse(start(x)),
                "\nEnd =", deparse(end(x)),
                "\nFrequency =", deparse(fr.x), "\n")
        else
            cat("Time Series:\nStart =", format(tsp(x)[1L]),
                "\nEnd =", format(tsp(x)[2L]),
                "\nFrequency =", deparse(fr.x), "\n")
    }
    print(.preformat.ts(x, calendar, ...), quote = FALSE, right = TRUE, ...)
    invisible(x)
}

## To be used in a  format.ts():
.preformat.ts <- function(x, calendar, ...)
{
    fr.x <- frequency(x)
    if(missing(calendar))
	calendar <- any(fr.x == c(4,12)) && length(start(x)) == 2L
    ## sanity check
    Tsp <- tsp(x)
    if(is.null(Tsp)) stop("series is corrupt, with no 'tsp' attribute")
    nn <- 1 + round((Tsp[2L] - Tsp[1L]) * Tsp[3L])
    if(NROW(x) != nn) {
        warning(gettextf("series is corrupt: length %d with 'tsp' implying %d",
                         NROW(x), nn), domain=NA, call.=FALSE)
        calendar <- FALSE
    }
    if(NCOL(x) == 1) { # could be 1-col matrix
        if(calendar) {
            if(fr.x > 1) {
                dn2 <-
                    if(fr.x == 12) month.abb
                    else if(fr.x == 4) {
                        c("Qtr1", "Qtr2", "Qtr3", "Qtr4")
                    } else paste0("p", 1L:fr.x)
                if(NROW(x) <= fr.x && start(x)[1L] == end(x)[1L]) {
                    ## not more than one period
                    dn1 <- start(x)[1L]
                    dn2 <- dn2[1 + (start(x)[2L] - 2 + seq_along(x))%%fr.x]
                    x <- matrix(format(x, ...), nrow = 1L , byrow = TRUE,
                                dimnames = list(dn1, dn2))
                } else { # more than one period
                    start.pad <- start(x)[2L] - 1
                    end.pad <- fr.x - end(x)[2L]
                    dn1 <- start(x)[1L]:end(x)[1L]
                    x <- matrix(c(rep.int("", start.pad), format(x, ...),
                                  rep.int("", end.pad)), ncol =  fr.x,
                                byrow = TRUE, dimnames = list(dn1, dn2))
                }
            } else { ## fr.x == 1
                tx <- time(x)
                attributes(x) <- NULL
                names(x) <- tx
            }
        } else { ##-- no 'calendar' --
            attr(x, "class") <- attr(x, "tsp") <- attr(x, "na.action") <- NULL
        }
    } else { # multi-column matrix
        rownames(x) <-
	    if(calendar && fr.x > 1) {
		tm <- time(x)
		t2 <- 1 + round(fr.x*((tm+0.001) %%1))
                ## protect people against themselves if they set options(digits=2)
                p1 <- format(floor(zapsmall(tm, digits = 7))) # yr
		if(fr.x == 12)
		    paste(month.abb[t2], p1)
		else
		    paste(p1, if(fr.x == 4) c("Q1", "Q2", "Q3", "Q4")[t2]
			  else format(t2))
	    } else
		format(time(x))
        attr(x, "class") <- attr(x, "tsp") <- attr(x, "na.action") <- NULL
    }
    x
}## {.preformat.ts}

plot.ts <-
    function (x, y = NULL, plot.type = c("multiple", "single"),
	      xy.labels, xy.lines, panel = lines, nc, yax.flip = FALSE,
	      mar.multi = c(0, 5.1, 0, if(yax.flip) 5.1 else 2.1),
	      oma.multi = c(6, 0, 5, 0), axes = TRUE, ...)
{
    plotts <-
	function (x, y = NULL, plot.type = c("multiple", "single"),
		  xy.labels, xy.lines, panel = lines, nc,
		  xlabel, ylabel, type = "l", xlim = NULL, ylim = NULL,
		  xlab = "Time", ylab, log = "", col = par("col"), bg = NA,
		  pch = par("pch"), cex = par("cex"),
		  lty = par("lty"), lwd = par("lwd"),
		  axes = TRUE, frame.plot = axes, ann = par("ann"),
                  cex.lab = par("cex.lab"), col.lab = par("col.lab"),
                  font.lab = par("font.lab"), cex.axis = par("cex.axis"),
                  col.axis = par("col.axis"), font.axis = par("font.axis"),
		  main = NULL, ...)
    {
	plot.type <- match.arg(plot.type)
	nser <- NCOL(x)

	if(plot.type == "multiple" && nser > 1) {
	    addmain <- function(main, cex.main=par("cex.main"),
				font.main=par("font.main"),
				col.main=par("col.main"), ...)
		## pass 'cex.main' etc	via "..." from main function
		mtext(main, side=3, line=3,
		      cex=cex.main, font=font.main, col=col.main, ...)
	    panel <- match.fun(panel)
	    nser <- NCOL(x)
	    if(nser > 10) stop("cannot plot more than 10 series as \"multiple\"")
	    if(is.null(main)) main <- xlabel
	    nm <- colnames(x) %||% paste("Series", 1L:nser)
	    if(missing(nc)) nc <- if(nser > 4) 2 else 1
	    nr <- ceiling(nser/nc)

	    oldpar <- par(mar = mar.multi, oma = oma.multi, mfcol = c(nr, nc))
	    on.exit(par(oldpar))
	    for(i in 1L:nser) {
		plot.default(x[, i], axes = FALSE, xlab="", ylab="",
                             log = log, col = col, bg = bg, pch = pch, ann = ann,
                             type = "n", ...)
		panel(x[, i], col = col, bg = bg, pch = pch,
		      cex = cex, lwd = lwd, lty = lty, type = type, ...)
		if(frame.plot) box(...)
		y.side <- if (i %% 2 || !yax.flip) 2 else 4
		do.xax <- i %% nr == 0 || i == nser
		if(axes) {
		    axis(y.side, xpd = NA, cex.axis = cex.axis,
			 col.axis = col.axis, font.axis = font.axis, ...)
		    if(do.xax)
			axis(1, xpd = NA, cex.axis = cex.axis,
			     col.axis = col.axis, font.axis = font.axis, ...)
		}
		if(ann) {
		    mtext(nm[i], y.side, line=3, cex=cex.lab, col=col.lab,
                          font=font.lab, ...)
		    if(do.xax)
			mtext(xlab, side=1, line=3, cex=cex.lab, col=col.lab,
			      font=font.lab, ...)
		}
	    }
	    if(ann && !is.null(main)) {
		par(mfcol=c(1,1))
		addmain(main, ...)
	    }
	    return(invisible())
	}
	## end of multiple plot section

	x <- as.ts(x)
	if(!is.null(y)) {
	    ## want ("scatter") plot of y ~ x
	    y <- hasTsp(y)
	    if(NCOL(x) > 1 || NCOL(y) > 1)
		stop("scatter plots only for univariate time series")
	    if (is.ts(x) && is.ts(y)) {
		xy <- ts.intersect(x, y)
		xy <- xy.coords(xy[,1], xy[,2], xlabel, ylabel, log)
	    } else
		xy <- xy.coords(x, y, xlabel, ylabel, log)
	    if(missing(xlab)) xlab <- xy$xlab
	    if(missing(ylab)) ylab <- xy$ylab
	    if(is.null(xlim)) xlim <- range(xy$x[is.finite(xy$x)])
	    if(is.null(ylim)) ylim <- range(xy$y[is.finite(xy$y)])
	    n <- length(xy $ x)
	    if(missing(xy.labels)) xy.labels <- (n <= 150)
	    do.lab <-
		if(is.logical(xy.labels))
		    xy.labels
		else {
		    if(!is.character(xy.labels))
			stop("'xy.labels' must be logical or character")
		    TRUE
		}
	    ptype <- if(do.lab) "n" else if(missing(type)) "p" else type

            dev.hold(); on.exit(dev.flush())
	    plot.default(xy, type = ptype,
			 xlab = xlab, ylab = ylab,
			 xlim = xlim, ylim = ylim, log = log, col = col, bg = bg,
			 pch=pch, cex=cex, lty=lty, lwd=lwd,
                         axes = axes, frame.plot = frame.plot,
			 ann = ann, main = main, ...)
	    if(missing(xy.lines)) xy.lines <- do.lab
	    if(do.lab)
		text(xy, labels =
		     if(is.character(xy.labels)) xy.labels
		     else if(all(tsp(x) == tsp(y)))
                         formatC(unclass(time(x)), width = 1)
		     else seq_along(xy$x),
		     col = col, cex = cex)
	    if(xy.lines)
		lines(xy, col = col, lty = lty, lwd = lwd,
		      type = if(do.lab) "c" else "l")
	    return(invisible())
	}
	## Else : no y, only x

	if(missing(ylab)) {
	    ylab <- colnames(x)
	    if(length(ylab) != 1L)
		ylab <- xlabel
	}
	## using xy.coords() mainly for the log treatment
	if(is.matrix(x)) {
	    k <- ncol(x)
	    tx <- time(x)
	    xy <- xy.coords(x = matrix(rep.int(tx, k), ncol = k),
			    y = x, log = log, setLab = FALSE)
	    xy$x <- tx
	}
	else xy <- xy.coords(x, NULL, log = log, setLab = FALSE)
	if(is.null(xlim)) xlim <- range(xy$x)
	if(is.null(ylim)) ylim <- range(xy$y[is.finite(xy$y)])
	plot.new()
	plot.window(xlim, ylim, log, ...)
	if(is.matrix(x)) {
	    for(i in seq_len(k))
		lines.default(xy$x, x[,i],
			      col = col[(i-1L) %% length(col) + 1L],
			      lty = lty[(i-1L) %% length(lty) + 1L],
			      lwd = lwd[(i-1L) %% length(lwd) + 1L],
			      bg  = bg [(i-1L) %% length(bg) + 1L],
			      pch = pch[(i-1L) %% length(pch) + 1L],
			      cex = cex[(i-1L) %% length(cex) + 1L],
			      type = type)
	}
	else {
	    lines.default(xy$x, x, col = col[1L], bg = bg, lty = lty[1L],
			  lwd = lwd[1L], pch = pch[1L],
			  cex = cex[1L], type = type)
	}
	if (ann)
	    title(main = main, xlab = xlab, ylab = ylab, ...)
	if (axes) {
	    axis(1, ...)
	    axis(2, ...)
	}
	if (frame.plot) box(...)
    }## {plotts}

    xlabel <- if (!missing(x)) deparse1(substitute(x))# else NULL
    ylabel <- if (!missing(y)) deparse1(substitute(y))
    plotts(x = x, y = y, plot.type = plot.type,
	   xy.labels = xy.labels, xy.lines = xy.lines,
	   panel = panel, nc = nc, xlabel = xlabel, ylabel = ylabel,
           axes = axes, ...)
}

lines.ts <- function(x, ...)
    lines.default(time(as.ts(x)), x, ...)


window.default <- function(x, start = NULL, end = NULL,
                           frequency = NULL, deltat = NULL,
                           extend = FALSE, ts.eps = getOption("ts.eps"), ...)
{
    x <- hasTsp(x)
    xtsp <- tsp(x)
    xfreq <- xtsp[3L]
    xtime <- time(x)

    if(!is.null(frequency) && !is.null(deltat) &&
       abs(frequency*deltat - 1) > ts.eps)
        stop("'frequency' and 'deltat' are both not NULL and are inconsistent")
    noFreq <- is.null(frequency) && is.null(deltat)
    yfreq <- if(noFreq)
                 xfreq
             else if(is.null(deltat))
                 frequency
             else if(is.null(frequency))
                 1/deltat
    ## fuzz to use for (start, end) boundary comparisons:
    eps_ <- ts.eps/xfreq
    thin <- round(xfreq/yfreq)
    if (yfreq > 0 && abs(xfreq/yfreq -thin) < ts.eps) {
        yfreq <- xfreq/thin
    } else if(!noFreq) {
        thin <- 1
        yfreq <- xfreq
            warning("'frequency' not changed")
    }
    start <-
        if(is.null(start))
            xtsp[1L]
        else switch(length(start),
                    start,
                    start[1L] + (start[2L] - 1)/xfreq,
                    stop("bad value for 'start'"))
    if(start < xtsp[1L]-eps_ && !extend) {
	start <- xtsp[1L]
	warning("'start' value not changed")
    }

    end <-
        if(is.null(end))
            xtsp[2L]
        else switch(length(end),
                    end,
                    end[1L] + (end[2L] - 1)/xfreq,
                    stop("bad value for 'end'"))
    if(end > xtsp[2L]+eps_ && !extend) {
	end <- xtsp[2L]
	warning("'end' value not changed")
    }

    if(start > end + eps_)
	stop("'start' cannot be after 'end'")

    if(!extend) {
        if(all(abs(start - xtime) > eps_))
            start <- xtime[(xtime > start) & ((start + 1/xfreq) > xtime)]

        if(all(abs(end - xtime) > eps_))
            end <- xtime[(xtime < end) & ((end - 1/xfreq) < xtime)]

        i <- seq.int(trunc((start - xtsp[1L]) * xfreq + 1.5),
                     trunc((end   - xtsp[1L]) * xfreq + 1.5), by = thin)
        y <- if(is.matrix(x)) x[i, , drop = FALSE] else x[i]
        ystart <- xtime[i[1L]]
        yend <- xtime[i[length(i)]]
        attr(y, "tsp") <- c(ystart, yend, yfreq)
    } else {
        ## first adjust start and end to the time base
        ## try to ensure that they are exactly n/xfreq
        stoff <- ceiling((start - xtsp[1L]) * xfreq - ts.eps)
        enoff <-   floor((end   - xtsp[2L]) * xfreq + ts.eps)
        ystart <- stoff/xfreq + xtsp[1L]
        yend   <- enoff/xfreq + xtsp[2L]
        # Rounding can cause problems #PR13272
        if (ystart > yend && (ystart - yend)*xfreq < ts.eps)
            yend <- ystart
        nold <- round(xfreq*(xtsp[2L] - xtsp[1L])) + 1
        ## both start and end could be outside time base
        ## and indeed the new ad old ranges might not intersect.
        i <- if(start > xtsp[2L]+eps_ || end < xtsp[1L] - eps_)
            rep(nold+1, floor(1+(end-start)*xfreq + ts.eps))
        else {
            i0 <- 1+max(0, stoff); i1 <- nold + min(0, enoff)
            c(rep.int(nold+1, max(0, -stoff)),
              if(i0 <= i1) i0:i1,
              rep.int(nold+1, max(0, enoff)))
        }
        y <- if(is.matrix(x)) rbind(x, NA)[i, , drop = FALSE] else c(x, NA)[i]
        attr(y, "tsp") <- c(ystart, yend, xfreq)
        if(yfreq != xfreq) y <- Recall(y, frequency = yfreq)
    }
    y
}

window.ts <- function (x, ...) as.ts(window.default(x, ...))

`window<-` <- function(x, ..., value) UseMethod("window<-")

`window<-.ts` <- function(x, start, end, frequency, deltat, ..., value)
{
    xtsp <- tsp(x)
    m <- match.call(expand.dots = FALSE)
    m$value <- NULL
    m$extend <- TRUE
    m$x <- x # do not attempt to re-evaluate *tmp* in replacement call
    ## need stats:: for non-standard evaluation
    m[[1L]] <- quote(stats::window)
    xx <- eval.parent(m)
    xxtsp <- tsp(xx)
    start <- xxtsp[1L]; end <- xxtsp[2L]
    if(start > end) stop("'start' > 'end'")
    if (start < xtsp[1L] || end > xtsp[2L]) {
        warning("extending time series when replacing values", call. = FALSE)
        x <- window(x, min(start, xtsp[1L]), max(end, xtsp[2L]), extend = TRUE)
    }
    xfreq <- xtsp[3L]
    xtimes <- round(xfreq*time(x))
    xxtimes <- round(xfreq * time(xx))

    ind <- match(xxtimes, xtimes)
    if(anyNA(ind)) stop("times to be replaced do not match")

    len <- length(ind)
    val_len <- NROW(value)
    if(!val_len) stop("no replacement values supplied")
    if(val_len > len) stop("too many replacement values supplied")
    if(val_len > 1L && (len %% val_len))
        stop("number of values supplied is not a sub-multiple of the number of values to be replaced")
    if(NCOL(x) == 1L) x[ind] <- value else x[ind, ] <- value
    x
}

`[.ts` <- function (x, i, j, drop = TRUE) {
    y <- NextMethod("[")
    if (missing(i))
	ts(y, start = start(x), frequency = frequency(x))
    else y
}

`[<-.ts` <- function (x, i, j, value) {
    y <- NextMethod("[<-")
    if (NROW(y) != NROW(x)) stop("only replacement of elements is allowed")
    y
}

t.ts <- function(x) {
    cl <- oldClass(x)
    other <- !(cl %in% c("ts","mts"))
    class(x) <- if(any(other)) cl[other]
    attr(x, "tsp") <- NULL
    t(x)
}

ts.plot <- function(..., gpars = list())
{
    dots <- list(...)
    pars <- c("xlab", "ylab", "xlim", "ylim", "col", "lty", "lwd",
              "type", "main", "sub", "log")
    m <- names(dots) %in% pars
    if(length(m)) {
        gpars <- c(gpars, dots[m])
        dots <- dots[!m]
    }
    sers <- do.call("ts.union", dots)
    if(is.null(gpars$ylab))
        gpars$ylab <- if(NCOL(sers) > 1) "" else deparse1(substitute(...))
    do.call("plot.ts", c(list(sers, plot.type = "single"), gpars))
}

arima.sim <- function(model, n, rand.gen = rnorm,
                      innov = rand.gen(n, ...), n.start = NA,
                      start.innov = rand.gen(n.start, ...), ...)
{
    if(!is.list(model)) stop("'model' must be list")
    if(n <= 0L) stop("'n' must be strictly positive")
    p <- length(model$ar)
    if(p) {
        minroots <- min(Mod(polyroot(c(1, -model$ar))))
        if(minroots <= 1) stop("'ar' part of model is not stationary")
    }
    q <- length(model$ma)
    if(is.na(n.start)) n.start <- p + q +
        ifelse(p > 0, ceiling(6/log(minroots)), 0)
    if(n.start < p + q) stop("burn-in 'n.start' must be as long as 'ar + ma'")
    d <- 0
    if(!is.null(ord <- model$order)) {
        if(length(ord) != 3L) stop("'model$order' must be of length 3")
        if(p != ord[1L]) stop("inconsistent specification of 'ar' order")
        if(q != ord[3L]) stop("inconsistent specification of 'ma' order")
        d <- ord[2L]
        if(d != round(d) || d < 0)
            stop("number of differences must be a positive integer")
    }
    if(!missing(start.innov) && length(start.innov) < n.start)
        stop(sprintf(ngettext(n.start,
                              "'start.innov' is too short: need %d point",
                              "'start.innov' is too short: need %d points"),
                     n.start), domain = NA)
    x <- ts(c(start.innov[seq_len(n.start)], innov[1L:n]), start = 1 - n.start)
    if(length(model$ma)) {
        x <- filter(x, c(1, model$ma), sides = 1L)
        x[seq_along(model$ma)] <- 0 # rather than NA
    }
    if(length(model$ar)) x <- filter(x, model$ar, method = "recursive")
    if(n.start > 0) x <- x[-(seq_len(n.start))]
    if(d > 0) x <- diffinv(x, differences = d)
    as.ts(x)
}