simulate.ets <- function(object, nsim=length(object$x), seed=NULL, future=TRUE, bootstrap=FALSE, innov=NULL, ...)
{
if(is.null(innov))
{
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else
{
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
}
if(is.null(tsp(object$x)))
object$x <- ts(object$x,frequency=1,start=1)
if(future)
initstate <- object$state[length(object$x)+1,]
else # choose a random starting point
initstate <- object$state[sample(1:length(object$x),1),]
if(bootstrap)
e <- sample(object$residuals,nsim,replace=TRUE)
else if(is.null(innov))
e <- rnorm(nsim,0,sqrt(object$sigma))
else if(length(innov)==nsim)
e <- innov
else
stop("Length of innov must be equal to nsim")
if(object$components[1]=="M")
e <- pmax(-1,e)
tmp <- ts(.C("etssimulate",
as.double(initstate),
as.integer(object$m),
as.integer(switch(object$components[1],"A"=1,"M"=2)),
as.integer(switch(object$components[2],"N"=0,"A"=1,"M"=2)),
as.integer(switch(object$components[3],"N"=0,"A"=1,"M"=2)),
as.double(object$par["alpha"]),
as.double(ifelse(object$components[2]=="N",0,object$par["beta"])),
as.double(ifelse(object$components[3]=="N",0,object$par["gamma"])),
as.double(ifelse(object$components[4]=="FALSE",1,object$par["phi"])),
as.integer(nsim),
as.double(numeric(nsim)),
as.double(e),
PACKAGE="forecast")[[11]],frequency=object$m,start=tsp(object$x)[2]+1/tsp(object$x)[3])
if(is.na(tmp[1]))
stop("Problem with multiplicative damped trend")
if(!is.null(object$lambda))
tmp <- InvBoxCox(tmp,object$lambda)
return(tmp)
}
# Simulate ARIMA model starting with observed data x
# Some of this function is borrowed from the arima.sim() function in the stats package.
# Note that myarima.sim() does simulation conditional on the values of observed x, whereas
# arima.sim() is unconditional on any observed x.
myarima.sim <- function (model, n, x, e, ...)
{
start.innov <- residuals(model)
innov <- e
data <- x
# Remove initial NAs
first.nonmiss <- which(!is.na(x))[1]
if(first.nonmiss > 1)
{
tsp.x <- tsp(x)
start.x <- tsp.x[1] + (first.nonmiss-1)/tsp.x[3]
x <- window(x,start=start.x)
start.innov <- window(start.innov, start=start.x)
}
model$x <- x
n.start <- length(x)
x <- ts(c(start.innov, innov), start = 1 - n.start, frequency=model$seasonal.period)
flag.noadjust <- FALSE
if(is.null(tsp(data)))
data <- ts(data,frequency=1,start=1)
if (!is.list(model))
stop("'model' must be list")
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)
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 (length(model$ma))
{
#MA filtering
x <- stats::filter(x, c(1, model$ma), method="convolution", sides = 1L)
x[seq_along(model$ma)] <- 0
}
##AR "filtering"
len.ar <- length(model$ar)
if(length(model$ar) & (len.ar <= length(data)))
{
if((model$seasonal.difference != 0) && (d != 0))
{
diff.data <- diff(data, lag=1, differences = d)
diff.data <- diff(diff.data, lag=model$seasonal.period, differences = model$seasonal.difference)
}
else if((model$seasonal.difference != 0) && (d == 0))
{
diff.data <- diff(data, lag=model$seasonal.period, differences = model$seasonal.difference)
}
else if((model$seasonal.difference == 0) && (d != 0))
diff.data <- diff(data, lag=1, differences = d)
else
diff.data <- data
x.new.innovations <- x[(length(start.innov)+1):length(x)]
x.with.data <- c(diff.data, x.new.innovations)
for(i in (length(diff.data)+1):length(x.with.data))
{
lagged.x.values <- x.with.data[(i-len.ar):(i-1)]
ar.coefficients <- model$ar[length(model$ar):1]
sum.mutliplied.x <- sum(lagged.x.values * (ar.coefficients))
x.with.data[i] <- x.with.data[i]+sum.mutliplied.x
}
x.end <- x.with.data[(length(diff.data)+1):length(x.with.data)]
x <- ts(x.end, start = 1, frequency=model$seasonal.period)
flag.noadjust <- TRUE
}
else if (length(model$ar)) # but data too short
{
#AR filtering for all other cases where AR is used.
x <- stats::filter(x, model$ar, method = "recursive")
}
if((d == 0) && (model$seasonal.difference == 0) && (flag.noadjust==FALSE)) # Adjust to ensure end matches approximately
{
# Last 20 diffs
if(n.start >= 20)
xdiff <- (model$x - x[1:n.start])[n.start-(19:0)]
else
xdiff <- model$x - x[1:n.start]
# If all same sign, choose last
if(all(sign(xdiff)==1) | all(sign(xdiff)==-1))
xdiff <- xdiff[length(xdiff)]
else # choose mean.
xdiff <- mean(xdiff)
x <- x + xdiff
}
if ((n.start > 0) && (flag.noadjust==FALSE))
{
x <- x[-(1:n.start)]
}
##
#####
#Seasonal undifferencing, if there is no regular differencing
if((model$seasonal.difference > 0) && (d == 0))
{
i <- length(data)-model$seasonal.difference*model$seasonal.period+1
seasonal.xi <- data[i:length(data)]
length.s.xi <- length(seasonal.xi)
x <- diffinv(x, lag=model$seasonal.period, differences=model$seasonal.difference, xi=seasonal.xi)[-(1:length.s.xi)]
data.new <- data
}
else
{
data.new <- data
}
###End seasonal undifferencing
##Regular undifferencing, if there is no seasonal differencing
if (d > 0 && (model$seasonal.difference == 0))
{
x <- diffinv(x, differences = d,xi=data.new[length(data.new)-(d:1)+1])[-(1:d)]
}
########
#Code for Undifferencing for where the differencing is both Seasonal and Non-Seasonal (Non-Seasonal First)
#Regular first
if((d > 0) && (model$seasonal.difference > 0))
{
delta.four <- diff(data, lag=model$seasonal.period, differences = model$seasonal.difference)
regular.xi <- delta.four[(length(delta.four)-model$seasonal.difference):length(delta.four)]
x <- diffinv(x, differences = d, xi=regular.xi[length(regular.xi)-(d:1)+1])[-(1:d)]
}
#Then seasonal
if((model$seasonal.difference > 0) && (d > 0))
{
i <- length(data)-model$seasonal.difference*model$seasonal.period+1
seasonal.xi <- data[i:length(data)]
length.s.xi <- length(seasonal.xi)
x <- diffinv(x, lag=model$seasonal.period, differences=model$seasonal.difference, xi=seasonal.xi)
x <- x[-(1:length.s.xi)]
data.new <- data
}
########
x <- ts(x[1:n],frequency=frequency(data),start=tsp(data)[2]+1/tsp(data)[3])
return(x)
}
simulate.Arima <- function(object, nsim=length(object$x), seed=NULL, xreg=NULL, future=TRUE, bootstrap=FALSE, innov=NULL, lambda=object$lambda, ...)
{
#Error check:
if(object$arma[7] < 0)
{
stop("Value for seasonal difference is < 0. Must be >= 0")
}
else if((sum(object$arma[c(3,4,7)])>0) && (object$arma[5] < 2))
{
stop("Invalid value for seasonal period")
}
####
#Random Seed Code
if(is.null(innov))
{
if (!exists(".Random.seed", envir = .GlobalEnv))
runif(1)
if (is.null(seed))
RNGstate <- .Random.seed
else
{
R.seed <- .Random.seed
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
}
#############End Random seed code
#Check for seasonal ARMA components and set flag accordingly. This will be used later in myarima.sim()
flag.s.arma <- (sum(object$arma[c(3,4)])>0)
#Check for Seasonality in ARIMA model
if(sum(object$arma[c(3,4,7)])>0)
{
#return(simulateSeasonalArima(object, nsim=nsim, seed=seed, xreg=xreg, future=future, bootstrap=bootstrap, ...))
if(sum(object$model$phi) == 0)
{
ar <- NULL
}
else
{
ar <- as.double(object$model$phi)
}
if(sum(object$model$theta) == 0)
{
ma <- NULL
}
else
{
ma <- as.double(object$model$theta)
}
order <- c(length(ar),object$arma[6],length(ma))
if(future)
{
model <- list(order=order, ar=ar, ma=ma,sd=sqrt(object$sigma2),residuals=residuals(object), seasonal.difference=object$arma[7], seasonal.period=object$arma[5], flag.seasonal.arma=flag.s.arma, seasonal.order=object$arma[c(3,7,4)])
}
else
{
model <- list(order=order, ar=ar, ma=ma,sd=sqrt(object$sigma2),residuals=residuals(object))
}
flag.seasonal.diff <- (object$arma[7] > 0)
}
else
{
####Non-Seasonal ARIMA specific code: Set up the model
order <- object$arma[c(1, 6, 2)]
if(order[1]>0)
ar <- object$model$phi[1:order[1]]
else
ar <- NULL
if(order[3]>0)
ma <- object$model$theta[1:order[3]]
else
ma <- NULL
if(object$arma[2] != length(ma))
stop("MA length wrong")
else if(object$arma[1] != length(ar))
stop("AR length wrong")
if(future)
{
model <- list(order=object$arma[c(1, 6, 2)],ar=ar,ma=ma,sd=sqrt(object$sigma2),residuals=residuals(object), seasonal.difference=0, flag.seasonal.arma=flag.s.arma, seasonal.order=c(0,0,0), seasonal.period=1)
}
else
{
model <- list(order=object$arma[c(1, 6, 2)],ar=ar,ma=ma,sd=sqrt(object$sigma2),residuals=residuals(object))
}
flag.seasonal.diff <- FALSE
###End non-seasonal ARIMA specific code
}
x <- object$x <- getResponse(object)
if(is.null(tsp(x)))
x <- ts(x,frequency=1,start=1)
if(!is.null(lambda))
x <- BoxCox(x,lambda)
n <- length(x)
d <- order[2]
if(bootstrap)
e <- sample(model$residuals,nsim+d,replace=TRUE)
else if(is.null(innov))
e <- rnorm(nsim+d, 0, model$sd)
else if(length(innov)==nsim)
e <- innov
else
stop("Length of innov must be equal to nsim")
use.drift <- is.element("drift", names(object$coef))
usexreg <- (!is.null(xreg) | use.drift)
if (!is.null(xreg))
{
xreg <- as.matrix(xreg)
if(nrow(xreg) < nsim)
stop("Not enough rows in xreg")
else
xreg <- xreg[1:nsim,]
}
if (use.drift)
{
dft <- as.matrix(1:nsim) + n
xreg <- cbind(dft,xreg)
}
narma <- sum(object$arma[1L:4L])
if(length(object$coef) > narma)
{
if (names(object$coef)[narma + 1L] == "intercept")
{
xreg <- cbind(intercept = rep(1, nsim), xreg)
object$xreg <- cbind(intercept = rep(1, n), object$xreg)
}
if(!is.null(xreg))
{
xm <- if (narma == 0)
drop(as.matrix(xreg) %*% object$coef)
else
drop(as.matrix(xreg) %*% object$coef[-(1L:narma)])
oldxm <- if(narma == 0)
drop(as.matrix(object$xreg) %*% object$coef)
else
drop(as.matrix(object$xreg) %*% object$coef[-(1L:narma)])
}
}
else
{
xm <- oldxm <- 0
}
if(future)
{
sim <- myarima.sim(model,nsim,x-oldxm,e=e) + xm
}
else
{
if(flag.seasonal.diff)
{
zeros <- object$arma[5]*object$arma[7]
sim <- arima.sim(model,nsim,innov=e)
sim <- diffinv(sim, lag=object$arma[5], differences=object$arma[7])[-(1:zeros)]
sim <- sim + xm
}
else
{
sim <- arima.sim(model,nsim,innov=e) + xm
}
}
if(!is.null(lambda))
sim <- InvBoxCox(sim,lambda)
return(sim)
}
simulate.ar <- function(object, nsim=object$n.used, seed=NULL, future=TRUE, bootstrap=FALSE, innov=NULL, ...)
{
if(is.null(innov))
{
if (!exists(".Random.seed", envir = .GlobalEnv))
runif(1)
if (is.null(seed))
RNGstate <- .Random.seed
else
{
R.seed <- .Random.seed
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
}
if(future)
{
model <- list(ar=object$ar,sd=sqrt(object$var.pred),residuals=object$resid, seasonal.difference=0, seasonal.period=1, flag.seasonal.arma=FALSE)
}
else
{
model <- list(ar=object$ar,sd=sqrt(object$var.pred),residuals=object$resid)
}
x.mean <- object$x.mean
object$x <- getResponse(object)
object$x <- object$x - x.mean
if(bootstrap)
e <- sample(na.omit(model$residuals),nsim,replace=TRUE)
else if(is.null(innov))
e <- rnorm(nsim, 0, model$sd)
else if(length(innov)==nsim)
e <- innov
else
stop("Length of innov must be equal to nsim")
if(future)
return(myarima.sim(model,nsim,x=object$x,e=e) + x.mean)
else
return(arima.sim(model,nsim,innov=e) + x.mean)
}
simulate.fracdiff <- function(object, nsim=object$n, seed=NULL, future=TRUE, bootstrap=FALSE, innov=NULL, ...)
{
x <- object$x <- getResponse(object)
# Strip initial and final missing values
xx <- na.ends(x)
n <- length(xx)
# Remove mean
meanx <- mean(xx)
xx <- xx - meanx
y <- undo.na.ends(x,diffseries(xx, d = object$d))
fit <- Arima(y, order = c(length(object$ar), 0, length(object$ma)),
include.mean = FALSE, fixed = c(object$ar, -object$ma))
# Simulate ARMA
ysim <- simulate(fit,nsim,seed,future=future,bootstrap=bootstrap,innov=innov)
# Undo differencing
return(unfracdiff(xx,ysim,n,nsim,object$d))
# bin.c <- (-1)^(0:(n + nsim)) * choose(object$d, (0:(n + nsim)))
# b <- numeric(n)
# xsim <- LHS <- numeric(nsim)
# RHS <- cumsum(ysim)
# bs <- cumsum(bin.c[1:nsim])
# b <- bin.c[(1:n) + 1]
# xsim[1] <- RHS[1] <- ysim[1] - sum(b * rev(xx))
# for (k in 2:nsim)
# {
# b <- b + bin.c[(1:n) + k]
# RHS[k] <- RHS[k] - sum(b * rev(xx))
# LHS[k] <- sum(rev(xsim[1:(k - 1)]) * bs[2:k])
# xsim[k] <- RHS[k] - LHS[k]
# }
# tspx <- tsp(x)
# return(ts(xsim,frequency=tspx[3],start=tspx[2]+1/tspx[3]))
}