# 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 == ""
##
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 ...
##
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])
##
## 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)
}