R/simulate.R

In ttnsdcn/forecast-package: Forecasting functions for time series

simulate.ets <- function(object, nsim=length(object$x), seed=NULL, future=TRUE, bootstrap=FALSE,...)
{
  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,f=1,s=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
    e <- rnorm(nsim,0,sqrt(object$sigma))
  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]],f=object$m,s=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,f=1,s=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 <- 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) && (d==0) && (model$seasonal.difference == 0) && (model$seasonal.order[1] != 0) && (model$seasonal.order[3] != 0) && (len.ar <= length(data))) 
	{
		##AR when D=0, d=0, SAR>0 and SMA>0
		x.new.innovations <- x[(length(start.innov)+1):length(x)]
		x.with.data <- c(data, x.new.innovations)
		for(i in (length(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(data)+1):length(x.with.data)]
		x <- ts(x.end, start = 1, frequency=model$seasonal.period)
		flag.noadjust <- TRUE
		
	} 
	else if(length(model$ar) && (model$seasonal.order[1] != 0) && (model$seasonal.order[3] != 0) && (len.ar <= length(data))) 
	{
		##AR when D>0, d>0, SAR>0 and SMA>0
		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)
		}
		
		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)) 
	{
		#AR filtering for all other cases where AR is used.
    x <- 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],f=frequency(data),s=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, ...)
{
	#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 (!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()
	if(sum(object$arma[c(3,4)])>0) 
	{
		flag.s.arma <- TRUE
	} 
	else 
	{
		flag.s.arma <- FALSE	
	}
	#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))
		}
		
		if(object$arma[7] > 0) 
		{
			flag.seasonal.diff <- TRUE
		} 
		else 
		{
			flag.seasonal.diff <- FALSE
		}
		
	} 
	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
	}
	

	
    if (is.element("x", names(object))) 
        x <- object$x
    else 
        x <- object$x <- eval.parent(parse(text = object$series))
		
	if(is.null(tsp(x)))
		x <- ts(x,f=1,s=1)

    n <- length(x)
	d <- order[2]
	if(bootstrap)
		e <- sample(model$residuals,nsim+d,replace=TRUE)
	else
		e <- rnorm(nsim+d, 0, model$sd)
		
	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(xreg, dft)
    }
	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(object$lambda))
		sim <- InvBoxCox(sim,object$lambda)
	return(sim)
}

simulate.ar <- function(object, nsim=object$n.used, seed=NULL, future=TRUE, bootstrap=FALSE, ...)
{
    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
    if (!is.element("x", names(object))) 
        object$x <- eval.parent(parse(text = object$series))
	if(is.null(tsp(object$x)))
		object$x <- ts(object$x,f=1,s=1)
    object$x <- eval.parent(parse(text = object$series)) - x.mean
	if(bootstrap)
		e <- sample(model$residuals,nsim,replace=TRUE)
	else
		e <- rnorm(nsim, 0, model$sd)
	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, ...)
{
	if (is.element("x", names(object))) 
        x <- object$x
    else 
		x <- object$x <- eval.parent(parse(text = as.character(object$call)[2]))
	if(is.null(tsp(x)))
		x <- ts(x,f=1,s=1)
    
    # 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)
	# 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,f=tspx[3],s=tspx[2]+1/tspx[3]))
}