######################################################################## ### sarfima.forecast: Forecast a time series using the SARFIMA Model ### ######################################################################## ### compute coefficients of (1-z)^m0(1-z^s)^ms difflag1coeffs.gen <- function(m0, n) c(choose(m0,0:m0)*(-1)^(0:m0),rep(0,n))[1:n] difflagScoeffs.gen <- function(ms, ss, n) c(as.vector(matrix(c(difflag1coeffs.gen(ms, ms+1), rep(0, (ss - 1) * (ms + 1))), ss, ms + 1, byrow = T)), rep(0, n))[1:n] difflag1Scoeffs.gen <- function(m0, ms, ss, n) conv(difflag1coeffs.gen(m0, n), difflagScoeffs.gen(ms, ss, n))[1:n] ### PMSE of BLP sarfima.PMSE.gen <- function(m0, ms, d0, ds, p0vec, q0vec, psvec, qsvec, ss, maxh, sigma2) { m <- m0 + ms * ss pijseq <- pij.SARFIMA.gen(d0, ds, p0vec, q0vec, psvec, qsvec, ss, maxh) y.psiseq <- inv.coef.gen(pijseq) m0ms.psiseq <- inv.coef.gen(difflag1Scoeffs.gen(m0, ms, ss, maxh)) x.psiseq <- (conv(m0ms.psiseq, y.psiseq))[1:maxh] list(PMSE=cumsum(x.psiseq^2) * sigma2, SRPMSE= sqrt(cumsum(x.psiseq^2) * sigma2)) } ### MAIN PROGRAM: sarfima.forecast ############################################################## # Forecast a time series using the SARFIMA Model. # # USAGE: sarfima.forecast(z, data, maxh, diag.plot = T, futuredata = NA, detail = F, con.level = 0.95) # # REQUIRED ARGUMENTS: # z a list of the output from sarfima.css # data a univariate time serires or a vector # maxh the number of time periods to forecast beyond the end of the series. # # OPTIONAL ARGUMENTS: # diag.plot logical flag: if TRUE, forecasts and predictive confidence itervals will be plotted. # futuredata future data (not to be used for estimating parameters). # detail logical flag: if TRUE, data, forecasts, predictive intervals, and so on will be provided. # con.level confidence level for prediction under the normality assumption. Default is 0.95. # # VALUE # a list containing the following elements: # timeseries$data same as input. # timeseries$futuredata same as input. # timeseries$xnhin forecasts of data used to compute the CSS function. # timeseries$xnh forecasts of future data. # upperCL upper confidence limits of forecasts designed by confidence level (con.level). # lowerCL lower confidence limits of forecasts designed by confidence level (con.level). # SRPMSE the approximate Square Root of Prediction Mean Squared Error. # PMSE approximate Prediction Mean Squared Error. # pred.error data (or futuredata) minus xnhin (or xnh). # r.pred.error 100 * pred.error / data (or futuredata) ################################################################################################# sarfima.forecast <- function(z, data, maxh, diag.plot = T, futuredata = NA, detail = F, con.level = 0.95) { arima.order <- z$model[[1]]$order seasonal.order <- z$model[[2]]$order p0 <- arima.order[1] m0 <- arima.order[2] q0 <- arima.order[3] ps <- seasonal.order[1] ms <- seasonal.order[2] qs <- seasonal.order[3] ss <- z$model[[2]]$period d0.est <- z$model[[3]]$d0.est ds.est <- z$model[[3]]$ds.est d0 <- z$initial.value$d0ds[1] ds <- z$initial.value$d0ds[2] xreg <- z$xreg x <- as.vector(data) xdiff0 <- if(m0 > 0) diff(x, lag = 1, differences = m0) else x yxreg <- if(ms > 0) diff(xdiff0, lag = ss, differences = ms) else xdiff0 mu <- nullto0(xreg) y0 <- yxreg - mu n <- length(y0) ### estimates ### resd0 <- z$estimates$d0 resds <- z$estimates$ds resar <- z$estimates$ar resma <- z$estimates$ma ressar <- z$estimates$sar ressma <- z$estimates$sma ressigma2 <- z$estimates$sigma2 n.d0 <- nullto0(resd0, b = d0) n.ds <- nullto0(resds, b = ds) n.ar <- nullto0(resar) n.ma <- nullto0(resma) n.sar <- nullto0(ressar) n.sma <- nullto0(ressma) ### residual and interpolation for y0 ### residual <- resid.SARFIMA.gen(n.d0, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n, y0) y0in <- c(y0[1], y0[2:n] - residual[2:n]) ### generate predictor of y0, y0nh ### pijseq <- pij.SARFIMA.gen(n.d0, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n + maxh) y0nh <- vector(length = maxh) yrevseq <- rev(y0) for(h in 1:maxh) { y0nh[h] <- - t(pijseq[2:(n + h)]) %*% yrevseq yrevseq <- c(y0nh[h], yrevseq) } ### generate predictor of yxreg, yxregnh, yxregin ### yxregnh <- y0nh + mu yxregin <- y0in + mu ### generate predictor of x and PMSEs ### m <- m0 + ms * ss if(m > 0) { aj <- - difflag1Scoeffs.gen(m0, ms, ss, m + 1)[2:(m + 1)] xnh <- vector(length = maxh) xrevseq <- rev(x) for(h in 1:maxh) { xnh[h] <- yxregnh[h] + t(aj) %*% xrevseq[1:m] xrevseq <- c(xnh[h], xrevseq) } xnhin <- vector(length = (m + n)) xnhin[1:m] <- x[1:m] xnhin[(m + 1):(m + n)] <- x[(m + 1):(m + n)] - residual } else { xnh <- yxregnh xnhin <- yxregin } y.psiseq <- inv.coef.gen(pijseq[1:maxh]) m0ms.psiseq <- inv.coef.gen(difflag1Scoeffs.gen(m0, ms, ss, maxh)) x.psiseq <- (conv(m0ms.psiseq, y.psiseq))[1:maxh] PMSEhseq <- cumsum(x.psiseq^2) * ressigma2 upperCL <- c(rep(NA, n + m), c(xnh + qnorm(1 - (1 - con.level) * 0.5) * sqrt(c(PMSEhseq)))) lowerCL <- c(rep(NA, n + m), c(xnh + qnorm((1 - con.level) * 0.5) * sqrt(c(PMSEhseq)))) if(is.rts(data)) { start.data <- start(data) freq.data <- frequency(data) end.data <- end(data) delta.data <- deltat(data) units.data <- units(data) } else { warning("Warning: data is not a regularly spaced time series\n") start.data <- 1 end.data <- n + m freq.data <- 1 delta.data <- 1 units.data <- NULL data <- rts(data, start = start.data, frequency = freq.data, units=units.data) } future <- rts(c(rep(NA, n + m - 1), as.vector(data[n + m]), as.vector(futuredata)), start = start.data, frequency = freq.data, units=units.data) xnhin <- rts(xnhin, start = start.data, frequency = freq.data, units=units.data) xnh <- rts(c(rep(NA, n + m), xnh), start = start.data, frequency = freq.data, units=units.data) upperCL <- rts(upperCL, start = start.data, frequency = freq.data, units=units.data) lowerCL <- rts(lowerCL, start = start.data, frequency = freq.data, units=units.data) SRPMSE <- rts(c(rep(NA, n + m), sqrt(PMSEhseq)), start = start.data, frequency = freq.data, units=units.data) ### tsplot ### if(diag.plot == T) { par(mfrow = c(2, 1)) zm <- max(n + m - max(2 * ss, 2 * maxh), floor((5 * (n + m))/6)) ts.plot(data, xnhin, xnh, upperCL, lowerCL, future, lty = c(1, 2, 3, 4, 4, 1), col = c(1, 6, 6, 8, 8, 1)) title("Plots of Data, Forecasts, and Forecast Intervals") ts.plot(data[ - c(1:zm)], xnhin[ - c(1:zm)], xnh[ - c(1:zm)], upperCL[ - c(1:zm)], lowerCL[ - c(1:zm)], future[ - c(1:zm)], lty = c(1, 2, 3, 4, 4, 1), col = c(1, 6, 6, 8, 8, 1)) title("Enlarged Plots of above") } ### output ### if(detail) { PMSE <- rts(c(rep(NA, n + m), PMSEhseq), start = start.data, frequency = freq.data, units=units.data) ### forecast error when futuredata is assigned ### pelength <- min(maxh, length(futuredata)) pred.error <- as.vector(futuredata)[1:pelength] - as.vector(xnh)[(n + m + 1):(n + m + pelength)] r.pred.error <- c(residual/as.vector(data)[(m + 1):(m + n)], pred.error/as.vector(futuredata)[1:pelength]) * 100 pred.error <- rts(c(rep(NA, m), residual, pred.error), start = start.data, frequency = freq.data, units=units.data) r.pred.error <- rts(c(rep(NA, m), r.pred.error), start = start.data, frequency = freq.data, units=units.data) futuredata <- rts(c(rep(NA, n + m), as.vector(futuredata)), start = start.data, frequency = freq.data, units=units.data) detail.list <- list(timeseries = ts.union(data, futuredata, xnhin, xnh, upperCL, lowerCL, PMSE, SRPMSE, pred.error, r.pred.error), datalist = list(start = start.data, end = end.data, freq = freq.data, deltat = delta.data), maxh = maxh, con.level = con.level) } else detail.list <- NULL list(forecast = xnh[(n + m + 1):(n + m + maxh)], SRPMSE = SRPMSE[(n + m + 1):(n + m + maxh)], detail = detail.list) } ### example 1 ################################################################################ data <- rts(arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 200), start = c(1980,1), frequency = 12) model <- list(list(order=c(1,0,1)), list(order=c(0,0,0), period=12), list(d0.est=T, ds.est=F)) fit1 <- sarfima.css(data, model) fore1 <- sarfima.forecast(fit1, data, maxh=24) ############################################################################################## ### example 2 ################################################################################ data <- rts(arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 264), start = c(1980,1), frequency = 12) est.data <- data[1:240] future.data <- data[241:264] model <- list(list(order=c(1,0,1)), list(order=c(0,0,0), period=12), list(d0.est=T, ds.est=F)) fit2 <- sarfima.css(est.data, model) fore2 <- sarfima.forecast(fit2, est.data, maxh=24, futuredata = future.data) ##############################################################################################