Nothing
## Do not edit this file manually.
## It has been automatically generated from *.org sources.
coreXarmaFilter <- function(x, eps, ar = numeric(0), ma = numeric(0),
p = length(ar), q = length(ma), n = length(x),
from = max(p, q) + 1, intercept = 0){
ceps <- if(q > 0){
filter(eps, c(1, ma), sides = 1L) # prepend '1' to the coef's.
# 2016-11-20 commenting out since:
# (a) filter() sets ceps[1:q] to NA (but this means that
# if p > 0 and p < q and from < q the result will be NA's,
# which is likely to be a programming error).
# (b) if from > q, ceps[1:q] is ignored below, so these values don't
# matter in this case.
# ceps[1:q] <- 0
}else{
eps
}
# intercept can be a constant or a vector of same length as ceps (and x)
ceps <- ceps + intercept
if(p > 0){
## this assumes that from > p
init <- x[(from - 1):(from - p)] # filter() requires init in reverse order;
# for "recursive", don't prepend '1' to 'ar'
ceps[from:n] <- filter(ceps[from:n], ar, method = "recursive", init = init)
}
x[from:n] <- ceps[from:n]
x
}
xarmaFilter <- function(model, x = NULL, eps = NULL, from = NULL, whiten = FALSE,
xcenter = NULL, xintercept = NULL){
# 2014-02-01 changing to "[[" to avoid partial matching for "p", etc.
# which was a long standing cause of puzzling errors.
ar <- model[["ar", exact = TRUE]] # model$ar
ma <- model[["ma", exact = TRUE]] # model$ma
center <- model[["center", exact = TRUE]] # model$center
intercept <- model[["intercept", exact = TRUE]] # model$intercept
# TODO: replace the above with:
# list2env(model, envir = environment()) # but need to change some is.null() below
# (no, it may be dangerous in a user facing function.)
p <- length(ar)
q <- length(ma)
if(is.null(from)) # set a value for "from" if not supplied
from <- 1 + max(p,q)
else if(from <= max(p,q))
stop("'from' must be greater than max(p,q)")
ct <- if(is.null(intercept)) 0 else intercept
if(!is.null(xintercept))
ct <- xintercept + ct
if(is.null(eps)){ # set initial values of eps to zero. copy 'x' to keep the attributes
# of x, eg if x is a "ts" object, and ensure consistency of returned
# value in the two cases for 'whiten'
eps <- x
eps[] <- 0 # makes sense mostly when whiten = TRUE
}else if(is.null(x)){
x <- eps
x[] <- 0
}
if(is.null(center) && is.null(xcenter)){
flag.center <- FALSE
y <- x
}else{
flag.center <- TRUE
mu <- if(is.null(xcenter)) 0 else xcenter
if(!is.null(center))
mu <- mu + center
y <- x - mu
}
n <- length(x)
stopifnot(n > 0)
if(whiten){ # compute residuals (whiten the series x if x is arma with these params)
eps <- coreXarmaFilter(x = eps, eps = y, ar = - ma, ma = - ar,
p = q, q = p, n = n, from = from, intercept = - ct)
eps
}else{ # compute a series from given residuals ("colour"/unwhiten eps)
y <- coreXarmaFilter(x = y, eps = eps, ar = ar, ma = ma,
p = p, q = q, n = n, from = from, intercept = ct )
x[from:n] <- if(flag.center) (y+mu)[from:n] # add back the center
else y[from:n]
x
}
}
# 2017-05-26: (somewhat) incompatible change - now uses prepareSimSarima()
sim_sarima <- function(model, n = NA, rand.gen = rnorm,
n.start = NA, x, eps, xcenter = NULL, xintercept = NULL, ...){
## decide if this patch should be permanent.
if(is.list(model) && is.null(model$sigma2))
model$sigma2 <- 1
fu <- prepareSimSarima(model = model, x = x, eps = eps, n = n,
n.start = n.start, xintercept = xintercept,
rand.gen = rand.gen)
res <- fu()
## TODO: 2020-03-02 sim_sarima is documented to return an object from class "ts"
## but doesn't!
## eg:
## frequency <- (environment(fu))$model$nseasons
## start <- if(is.null(frequency) || frequency == 1) 1 else frequency
## res <- ts(res, start = start, frequency = frequency)
##
## But it may be better to add arguments for the class, start, end, ...
res
}
model2filter <- function(model){
if(!is(model, "VirtualFilterModel"))
model <- do.call("new", c(list("SarimaModel"), model))
co <- modelCoef(model, "ArmaModel")
n_ur <- nUnitRoots(model)
list(ar = co$ar, ma = co$ma, # get 'ar' and 'ma' from 'co'
sigma2 = model@sigma2, # the rest from 'model'; could use sigmaSq(model), etc.
center = model@center,
intercept = model@intercept,
n_ur = n_ur # 2018-07-28 new component
)
}
## helper function - ensure that the argument is a list with components
## before, init and main
.beforeInitMain <- function(x){
res <- list(before = numeric(0), init = numeric(0), main = numeric(0))
if(missing(x) || is.null(x))
res
else if(inherits(x, "list")){
# x[c("before", "init", "main")] would return elements named <NA>
# for absent components
if(!is.null(x$before)) res$before <- x$before
if(!is.null(x$init)) res$init <- x$init
if(!is.null(x$main)) res$main <- x$main
res
}else{
res$main <- x
res
}
}
prepareSimSarima <- local({
## these variables will be shadowed by assignments in the body of prepareSimSarima
## they are defined for somewhat more clear code; 'R CMD check' is also happier.
x <- eps <- NULL
ar <- ma <- numeric(0)
p <- q <- 0
from <- NULL # TODO: rethink!
# the default for n.start is NA, so that in future automatic choice may be offered
function(model, x = NULL, eps = NULL, n, n.start = NA, xintercept = NULL,
rand.gen = rnorm){
eps <- .beforeInitMain(eps)
x <- .beforeInitMain(x)
ct <- .beforeInitMain(xintercept)
if(length(eps$before) != length(x$before)){
if(length(eps$before) == 0)
eps$before <- numeric(length(x$before))
else if(length(x$before) == 0)
x$before <- numeric(length(eps$before))
else
stop("Lengths of xbefore and innovbefore must be equal if both are present.")
}
if(length(eps$init) != length(x$init)){
if(length(eps$init) == 0)
eps$init <- numeric(length(x$init))
else if(length(x$init) == 0)
x$init <- numeric(length(eps$init))
else
stop("Lengths of xinit and innovinit must be equal if both are present.")
}
n.before <- length(x$before)
if(is.na(n.start))
n.start <- 0
if(missing(n) || is.null(n))
n <- length(eps$main) + length(eps$init) - n.start
else{ # adjust x$main and eps$main
eps$main <- c(eps$main,
numeric(n.start + n - length(eps$init) - length(eps$main)) )
x$main <- c(x$main,
numeric(n.start + n - length(x$init) - length(x$main)) )
}
from <- 1 + length(x$before) + length(x$init)
n.drop <- n.before + n.start
## TODO: use new variables further below; for now save the old values in eps0,x0
eps0 <- eps
x0 <- x
eps <- c(eps$before, eps$init, eps$main)
x <- c(x$before, x$init, numeric(n.start + n - length(x$init)))
## instead, could use:
## x <- c(x$before, x$init, x$main)
## but it is preferable to state explicitly that any values in x$main
## are discarded.
## this is guaranteed by the above code but state the intention:
stopifnot(length(eps) == length(x))
ctt <- if(length(ct$main) == 0)
0
else if(length(ct$main) == 1){ # constant xintercept
ct$main
}else if(length(ct$main) > 1){ # non-constant xintercept
if(length(ct$before) == 0) # TODO: NA on place of 0?
ct$before <- numeric(length(eps0$before))
if(length(ct$init) == 0)
ct$init <- numeric(length(eps0$init))
c(ct$before, ct$init, ct$main)
}
stopifnot(length(x) == length(eps),
length(ctt) == 1 || length(ctt) == length(eps)
)
filtmodel <- model2filter(model)
ar <- filtmodel$ar
ma <- filtmodel$ma
mu <- filtmodel$center
intercept <- filtmodel$intercept
sigma <- sqrt(filtmodel$sigma2)
p <- length(ar)
q <- length(ma)
if(intercept != 0)
ctt <- ctt + intercept
flag.center <- mu != 0
n0 <- n
rand.gen0 <- match.fun(rand.gen)
## 2018-07-28 krapka - for the case when no unit roots and no initial values
## => generate initial values
flag.stat.init.values <- filtmodel$n_ur == 0 && from == 1 # 2018-07-28
mpq <- max(p,q) ## TODO: what if mpq = 0?
if(flag.stat.init.values){
acv <- autocovariances(filtmodel[c("ar", "ma", "sigma2")], maxlag = mpq)
partial_sds <- sqrt(.comboAcvf(acv, "psigma2"))[-1]
partial_coefs <- acf2AR(acv[])
}
## TODO: tryabva po-palna proverka na razmernostite!
## TODO: drugi varianti, kontrolirani s argument or otherwise;
## e.g., use separate functions instead of conditionals
f <- function(n = n0, rand.gen, ...){
if(missing(rand.gen))
rand.gen <- rand.gen0
nx <- length(x)
## simulate (standardised) innovations and multiply by sigma
## 2022-03-24 was: eps.main <- sigma * rand.gen(n.start + n, ...)
eps <- if(from > 1) # 2018-07-28 added if()
c(eps[1:(from-1)], sigma * rand.gen(length(eps) - (from - 1), ...))
else
sigma * rand.gen(n.start + n, ...)
## TODO: may need to modify x if n != n0
if(flag.stat.init.values && mpq > 0){ # TODO: needs more thought
x[1:mpq] <- partial_sds * eps[1:mpq]
if(mpq > 1){
for(i in 2:mpq){
x[i] <- x[i] + sum(x[1:(i-1)] * partial_coefs[i-1, 1:(i-1)])
}
}
x[1:mpq] <- x[1:mpq] + mu # since mu will be subtracted below
}
y <- if(flag.center)
x - mu
else
x
stopifnot(length(x) == length(eps), mpq < length(x)) # 2022-03-24
y <- coreXarmaFilter(x = y, eps = eps, ar = ar, ma = ma, p = p, q = q, # n = n,
# from = from, TODO: temporary commenting out (v. 0.4-5)
intercept = ctt )
x[from:nx] <-
if(flag.center)
(y + mu)[from:nx]
else
y[from:nx]
#browser()
if(n.drop > 0) # 2018-07-28 added if()
x[-(1:n.drop)] # TODO: keep the class of x?
else
x # TODO: but this branch will keep it!
}
class(f) <- "simSarimaFun"
f
}
})
print.simSarimaFun <- function(x, ...){
print("A function for simulation of the following Sarima model.")
e <- environment(x)
print(e$model) # TODO: do this more friendly.
print("The parameters of the simulation are:")
print("TODO: this print method is unfinished.")
invisible(NULL)
}
sarima.f <- function(past = numeric(length(ar)),
n = max(2*length(past),12),
ar = numeric(0), ma = numeric(0),
intercept = 0,
pasteps = numeric(length(ma)), trend = numeric(n)){
p <- length(ar)
q <- length(ma)
m <- max(p,q)
if(length(past) > m) # keep only m past values
past <- past[length(past) + (-m+1):0]
res <- c(numeric(max(m-p,0)), past, numeric(n))
if(length(pasteps) > m) # keep only m past values
pasteps <- past[length(pasteps) + (-m+1):0]
eps <- c(numeric(max(m-q,0)), pasteps, numeric(n)) # more care needed
trend <- c(numeric(m), trend)
for(i in m + (1:n) ){
res[i] <- intercept + trend[i] + eps[i]
if(p>0)
for(j in 1:p){
res[i] <- res[i] + ar[j]*res[i-j]
}
if(q>0)
for(j in 1:q){
res[i] <- res[i] + ma[j]*eps[i-j]
}
}
res[-(1:m)] # Return the forecasts only.
}
fun.forecast <- function( past
, n=max(2*length(past),12)
#, trend = numeric(n)
, eps = numeric(n)
, pasteps
, ...
){
## model <- sarima.mod(...)
sarima <- new("SarimaModel", ...)
co <- filterCoef(sarima, "BD")
ar = co$ar
ma = co$ma
center = sarima@center
intercept = sarima@intercept
if(missing(past))
past <- numeric(length(ar))
if(missing(pasteps))
pasteps <- numeric(length(ma))
fullintercept <- intercept
if(center != 0 && sarima@iorder == 0 && sarima@siorder == 0) # Why only in this case?
# (ans: because otherwise (1 - sum(ar) ) = 0)
fullintercept <- intercept + (1-sum(ar)) * center
res <- sarima.f(past = past, n = n, ar = ar, ma = ma
, intercept = fullintercept
, pasteps = pasteps
) # trend=
ts(res) # Return the forecasts only. TODO: set period?
}
## ## for now use sim_sarima as default
## simSarima <- sim_sarima
##
## setGeneric("simSarima", signature = c("model"))
##
## 2016-10-26 commenting out temporarily after removal of arCoef and maCoef
##
## setMethod("simSarima", "ArmaModel",
## function(model, init = NULL, rand.gen = rnorm, info = "print", ...){
## mo <- sarima.mod(ar = arCoef(model), ma = maCoef(model), mean = mean(model))
## sim_sarima(model = mo, ...)
## }
## )
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.