Nothing
# These functions are handling the estimation of MAPA
# Nikolaos Kourentzes 2016
mapaest <- function(y, ppy=NULL, minimumAL=1, maximumAL=ppy, paral=c(0,1,2), display=c(0,1),
outplot=c(0,1), model="ZZZ", type=c("ets","es"), xreg=NULL,
pr.comp=0,...) {
# Estimate MAPA for a time series
#
# Inputs:
# y = In sample observations of a time series (vector)
# ppy = Periods in a season of the time series at the sampled frequency.
# If y is a ts object then this is taken from its frequency,
# unless overriden.
# minimumAL = Lowest aggregation level to use. Default = 1, maximumAL>1
# maximumAL = Highest aggregation level to use. Default = ppy
# paral = Use parallel processing. 0 = no; 1 = yes (requires initialised cluster);
# 2 = yes and initialise cluster. Default is 0.
# display = Display calculation progress in console. 0 = no; 1 = yes. Default is 0.
# outplot = Provide output plot. 0 = no; 1 = yes. Default is 1.
# model = Allow only that type of ETS at each aggregation level. This follows similar
# coding to the ets function. The first letter refers to the error
# type ("A", "M" or "Z"); the second letter refers to the trend
# type ("N","A","Ad","M","Md", "Z", "X" or "Y"); and the third letter
# refers to the season type ("N","A","M", "Z", "X" or "Y"). The letters mean:
# "N"=none, "A"=additive, "M"=multiplicative, "Z"=automatically selected,
# "X"=automatically select between none and additive and "Y"=automatically
# select between none and multiplicative. "X" and "Y" supported only by
# type=="es". If used with type=="ets" a warning will be given and they will
# default to "Z". A "d" for trend implies damped.
# By default model="ZZZ". If due to sample limitation ETS cannot be calculated
# at an aggregation level for the selected model, then no estimation is done
# for that specific level. For aggregation levels that seasonality becomes
# 1 then a non-seasonal model is estimated.
# type = What type of exponential smoothing implementation to use. "es" - use from
# the smooth package; "ets" use from the forecast package
# xreg = Vector or matrix of exogenous variables to be included in the MAPA.
# If matrix then rows are observations and columns are variables.
# Must be at least as long as in-sample. Additional observations are unused.
# Note that including xreg will force type="es".
# pr.comp = MAPAx can use principal component analysis to preprocess xreg. When comp is -1
# then the number of retained components is chosen automatically. When comp=0
# then no pre-processing is performed and the original xreg is used. Any other
# value represents the number of principal components retained.
# ... = Pass additional arguments to es and ets.
#
# Output:
# mapafit = Estimated MAPA model structure
# Default type
paral <- paral[1]
display <- display[1]
outplot <- outplot[1]
type <- type[1]
# Get ppy and maximumAL
if (is.null(ppy)){
if (isa(y,"ts")){
ppy <- frequency(y)
if (is.null(maximumAL)){maximumAL <- ppy}
} else {
stop(paste("Input ppy is not given and y input not ts class.",
"Please provide the periods in a season of the time series",
"at the sampled frequency."))
}
}
# If xreg is used then force type="es" and give a warning
if (!is.null(xreg) & type=="ets"){
type <- "es"
warning('Only type="es" can accept xreg inputs. Forcing appropriate type.')
}
# Convert xreg to array and check its size
if (!is.null(xreg)){
n <- length(y)
# Force xreg to be a matrix, each column a variable
if (!is.null(dim(xreg))){
p <- dim(xreg)[2]
} else {
p <- 1
}
xreg <- matrix(xreg,ncol=p)
# Check that it contains at least n observations and trim
xn <- dim(xreg)[1]
if (xn < n){
stop("Number of observations of xreg must be equal of exceed the length of y.")
} else {
xreg <- matrix(xreg[1:n,],ncol=p)
}
# Perform preprocessing
x.m <- colMeans(xreg)
x.sd <- apply(xreg,2,sd)
pr.comp <- list("pr.comp"=pr.comp, "mean"=x.m, "sd"=x.sd)
pr.out <- mapaprcomp(xreg,pr.comp)
xreg <- pr.out$x.out
pr.comp <- pr.out$pr.comp
} else {
# No xreg, so mark no pre-processing
pr.comp <- list("pr.comp"=0, "mean"=NULL, "sd"=NULL)
}
# Suggest to use es instead of ets when ppy > 24
if (ppy > 24 & type == "ets"){
warning('Your data have potential seasonality longer than 24 periods. Consider switching to type="es".')
}
# Make sure that maximumAL > 1 and larger than minimumAL
if (maximumAL == 1){
maximumAL = maximumAL + 1
}
if (minimumAL>=maximumAL){
stop("maximumAL must be larger than minimumAL")
}
# Setup parallel processing if required
if (paral == 2){
crs <- detectCores()
cl <- makeCluster(getOption("cl.cores", crs))
writeLines(paste("Running with", crs, 'cores'))
}
observations <- length(y) # number of observations for the in-sample data
# Override model if used with ETS
if (type=="ets"){
mn <- nchar(model)
model.set <- substring(model, seq(1,mn,1), seq(1,mn,1))
idx <- ((model.set)=="X" | (model.set)=="Y")
if (any(idx)){
warning('When type=="ets" neither "X" or "Y" can be used in model specification. These will be changed to "Z".')
model.set[idx] <- "Z"
model <- paste0(model.set,collapse="")
}
}
# Aggregation and estimation
if (paral != 0){ # Parallel run
mapafit <- clusterApplyLB(cl, 1:(maximumAL-minimumAL+1), mapaest.loop,
y=y, minimumAL=minimumAL, maximumAL=maximumAL, observations=observations, ppy=ppy,
display=display,model=model,type=type,xreg=xreg,pr.comp=pr.comp,...)
} else { # Serial run
mapafit <- vector("list", (maximumAL-minimumAL+1))
for (i in 1:(maximumAL-minimumAL+1)){
mapafit[[i]] <- mapaest.loop(i, y, minimumAL, maximumAL, observations,
ppy, display,model=model,type=type,xreg=xreg,pr.comp=pr.comp,...)
}
}
if (paral == 2){
# Stop parallel processing
stopCluster(cl)
}
# Process output
mapafit <- do.call(rbind, mapafit) # Re-arrange output for clusterApplyLB function
rownames(mapafit) <- paste("AL",minimumAL:maximumAL,sep="")
mapafit <- structure(mapafit, class = "mapa.fit")
# Plot model selection summary
if (outplot == 1){
plot(mapafit)
}
# Return output
return(mapafit)
}
#-------------------------------------------------
mapaest.loop <- function(ALi, y, minimumAL, maximumAL, observations,
ppy, display, model, type, xreg=NULL, pr.comp, ...){
# Internal function for running a single loop in mapaest
# Create ETS model strings
mn <- nchar(model)
model.set <- substring(model, seq(1,mn,1), seq(1,mn,1))
if (type == "ets"){
if (model.set[2]=="Z"){
ets.damped <- NULL
} else if (mn==4){
ets.damped <- TRUE
} else {
ets.damped <- FALSE
}
} else if (type == "es"){
if (mn == 4){
model.set <- c(model.set[1],paste0(model.set[2],model.set[3]),model.set[mn])
}
} else {
stop("Unrecognised exponential smoothing implementation type.")
}
ALvec <- minimumAL:maximumAL
AL <- ALvec[ALi]
# Console output
if (display==1){
txtc <- paste("Aggregation level: ",AL,"/",maximumAL,
" (",round(100*ALi/(maximumAL-minimumAL+1),2),"%)",sep="")
cat(txtc)
}
q <- observations %/% AL # observation in the aggregated level
# r <- observations %% AL # observation to discard from the beginning of the series
ppyA <- ppy %/% AL # periods per year for the aggregated level
if (ppy %% AL != 0){
ppyA <- 1
}
# Check if requested ETS model is possible for current AL
# Different treatment of damped trend and npars between packages
if (type == "ets"){
npars <- 2
if (!is.null(ets.damped)){
npars <- npars + as.numeric(ets.damped)}
} else if (type == "es"){
npars <- 3
if (nchar(model.set[2]) == 2){
npars <- npars + 3}
}
# Common amongst both types
if (model.set[2] == "A" | model.set[2] == "M"){
npars <- npars + 2}
if (tail(model.set,1) == "A" | tail(model.set,1) == "M"){
npars <- npars + ppyA}
if (q <= npars){
q <- 1} # This will not estimate current AL
# Adjust for explanatory variables
if (!is.null(xreg)){
xreg.p <- dim(xreg)[2]
} else {
xreg.p <- 0
}
if (q >= 4 + xreg.p){ # Fit only when there is enough sample
# Aggregation
yA <- colMeans(matrix(tail(y,q*AL),nrow=AL))
ats <- ts(yA, frequency = ppyA) # create the time series
# Old code
# yA <- array(0, dim=c(q)) # in-sample aggregated values will be saved here
# for (j in 1:q){ # calculate the aggregate values
# yA[j] <- mean(y[(r+1+(j-1)*AL):(r+j*AL)])
# }
# Aggregate xreg
if (!is.null(xreg)){
p <- dim(xreg)[2]
xregA <- array(NA, c(q,p))
for (k in 1:p){
xregA[,k] <- colMeans(matrix(tail(xreg[,k],q*AL),nrow=AL))
}
} else {
xregA <- NULL
}
# Check if seasonality exists and select appropriate model
if ((tail(model.set,1) == "A" | tail(model.set,1) == "M") && ppyA == 1){
mapa.model <- paste0(model.set[1],model.set[2],"N")
} else if ((ppyA == 2 && q <= 3*ppyA) & type == "ets"){
# This is added to deal with the bug in the fourier function that
# is used for the initialisation of ets
mapa.model <- paste0(model.set[1],model.set[2],"N")
} else {
mapa.model <- paste0(model.set[1],model.set[2],tail(model.set,1))
}
# Fit exponential smoothing
if (type == "ets"){
fit <- ets(ats,model=mapa.model,damped=ets.damped,...)
# If time series is constant then ets does not return the same
# structure, correct for that
if (length(names(fit)) != 18){
fit.temp = fit
fit <- list("loglik"=NULL,"aic"=NULL,"bic"=NULL,"aicc"=NULL,"mse"=NULL,
"amse"=NULL,"fit"=NULL,"residuals"=NULL,"fitted"=NULL,
"states"=NULL,"par"=NULL,"m"=NULL,"method"=NULL,
"components"=NULL,"call"=NULL,"initstate"=NULL,
"sigma2"=NULL,"x"=NULL)
fit.temp.names <- names(fit.temp)
fit.names <- names(fit)
for (fi in 1:length(fit.names)){
if (sum(fit.temp.names == fit.names[fi])>0){
eval(parse(text=paste0("fit$",fit.names[fi]," <- fit.temp$",fit.names[fi])))
}
}
}
} else if (type == "es"){
# When xreg is a single column array es wants the input as a vector
if (!is.null(xreg)){
if (dim(xregA)[2] ==1){
xregA <- as.vector(xregA)
}
}
# Turn off warnings for es - this is done when the model reduces pool
# due to sample size.
fittemp <- suppressWarnings(es(ats,model=mapa.model,silent=TRUE,xreg=xregA, ...))
acceptres <- c('model', 'timeElapsed', 'data', 'holdout', 'fitted',
'residuals', 'forecast', 'states', 'profile', 'profileInitial',
'persistence', 'phi', 'transition', 'measurement', 'initial',
'initialType', 'initialEstimated', 'orders', 'arma', 'constant',
'nParam', 'occurrence', 'formula', 'regressors', 'loss',
'lossValue', 'logLik', 'distribution', 'scale', 'other',
'B', 'lags', 'lagsAll', 'FI', 'ICs', 'call', 'bounds')
fit <- fittemp[unlist(lapply(acceptres,function(x){which(names(fittemp) == x)}))]
fit$smoothLength <- length(fit)
}
# Common part
fit$use <- TRUE
} else {
# No fit is taking place
if (type == "ets"){
fit <- list("loglik"=NULL,"aic"=NULL,"bic"=NULL,"aicc"=NULL,"mse"=NULL,
"amse"=NULL,"fit"=NULL,"residuals"=NULL,"fitted"=NULL,
"states"=NULL,"par"=NULL,"m"=NULL,"method"=NULL,
"components"=NULL,"call"=NULL,"initstate"=NULL,
"sigma2"=NULL,"x"=NULL,"use"=FALSE)
} else if (type == "es"){
fit <- list('model'=NULL, 'timeElapsed'=NULL, 'data'=NULL, 'holdout'=NULL,
'fitted'=NULL, 'residuals'=NULL, 'forecast'=NULL, 'states'=NULL,
'profile'=NULL, 'profileInitial'=NULL, 'persistence'=NULL,
'phi'=NULL, 'transition'=NULL, 'measurement'=NULL, 'initial'=NULL,
'initialType'=NULL, 'initialEstimated'=NULL, 'orders'=NULL,
'arma'=NULL, 'constant'=NULL, 'nParam'=NULL, 'occurrence'=NULL,
'formula'=NULL, 'regressors'=NULL, 'loss'=NULL, 'lossValue'=NULL,
'logLik'=NULL, 'distribution'=NULL, 'scale'=NULL, 'other'=NULL,
'B'=NULL, 'lags'=NULL, 'lagsAll'=NULL, 'FI'=NULL, 'ICs'=NULL,
'call'=NULL, 'bounds'=NULL,
'use'==FALSE)
}
}
fit$AL <- AL
fit$original.ppy <- ppy
if (type == "ets"){
fit$etstype <- "ets"
} else if (type == "es"){
fit$etstype <- "es"
}
fit$pr.comp <- pr.comp
# Update console display
if (display==1){
nc = nchar(txtc)
cat(rep("\r",nc))
cat(rep(" ",nc))
cat(rep("\r",nc))
}
# Return loop result
return(rbind(fit))
}
# #-------------------------------------------------
# plotmapa <- function(mapafit){
#
# warning("'plotmapa' is deprecated.\n Use plot() instead.")
# plot.mapa.fit(mapafit)
#
# }
#-------------------------------------------------
summary.mapa.fit <- function(object,...){
print(object)
}
#-------------------------------------------------
print.mapa.fit <- function(x,...){
# Get basic settings
idx.use <- which(colnames(x)=="use")
idx.AL <- which(colnames(x)=="AL")
idx.ppy <- which(colnames(x)=="original.ppy")
ets.type <- x[1,which(colnames(x)=="etstype")]
# Get settings from mapafit
ALs <- as.numeric(x[x[,idx.use]==TRUE, idx.AL])
minimumAL <- min(ALs)
maximumAL <- max(ALs)
ppy <- as.numeric(x[1,idx.ppy])
# Plot model selection summary
ALplot <- 1:(maximumAL-minimumAL+1)
ALplot <- ALplot[unlist(x[,idx.use])==TRUE]
cat(paste0("MAPA fitted using ",ets.type, "\tOriginal frequency: ",ppy,"\n"))
if (ets.type == "ets"){
mdls <- x[unlist(x[,idx.use]),13]
mdls <- lapply(mdls,function(x){gsub(',','',x)})
x.n <- rep(0,sum(unlist(x[,idx.use])))
} else if (ets.type == "es"){
mdls <- x[unlist(x[,idx.use]),1]
# Check for xreg
idx.x <- which(colnames(x)=="initial")
x.n <- lapply(x[,idx.x],function(x){length(x$xreg)})
x.n <- unlist(x.n[unlist(x[,idx.use])])
}
for (i in 1:(maximumAL-minimumAL+1)){
if (x.n[i] == 0){ # No xreg
cat(paste0("Aggregation level: ",ALplot[i],"\t","Method: ",mdls[[i]],"\n"))
}else {
cat(paste0("Aggregation level: ",ALplot[i],"\t","Method: ",mdls[[i]],"+X(",x.n[i],")","\n"))
}
}
}
#-------------------------------------------------
plot.mapa.fit <- function(x,xreg.plot=c(TRUE,FALSE),...){
# Produce estimated MAPA fit plot
#
# Inputs:
# mapafit = Fitted MAPA model (from mapaest)
# xreg.plot = Add infromation about xreg in the figure.
xreg.plot <- xreg.plot[1]
# Find locations in mapafit
idx.use <- which(colnames(x)=="use")
idx.AL <- which(colnames(x)=="AL")
idx.ppy <- which(colnames(x)=="original.ppy")
ets.type <- x[1,which(colnames(x)=="etstype")]
# Get settings from mapafit
ALs <- as.numeric(x[x[,idx.use]==TRUE, idx.AL])
minimumAL <- min(ALs)
maximumAL <- max(ALs)
ppy <- as.numeric(x[1,idx.ppy])
# Get permitted seasonal levels
perm_seas <- (ppy %% minimumAL:maximumAL) == 0 & (minimumAL:maximumAL < ppy)
# Plot model selection summary
ALplot <- 1:(maximumAL-minimumAL+1)
ALplot <- ALplot[unlist(x[,idx.use])==TRUE]
layout(matrix(1, 1, 1, byrow = TRUE))
comps <- array(0,c(max(ALplot),4))
comps.char <- array(0,c(max(ALplot),3))
x.n <- vector("numeric",max(ALplot))
for (AL in 1:max(ALplot)){
# Get fitted components at each AL
if (ets.type == "ets"){
components <- x[[AL, 14]]
if (components[[4]] == "TRUE"){
components[[4]] <- "d"
} else {
components[[4]] <- ""
}
components <- components[c(1,2,4,3)]
components <- c(components[1],paste0(components[2],components[3]),tail(components,1))
ttl <- "ETS components"
# Check for xreg
x.n[AL] <- 0
} else if (ets.type == "es"){
# Get ES components
components <- x[[AL, 1]]
components <- strsplit(components,"")
locstart <- which(unlist(lapply(components,function(x){x=="("})))
locend <- which(unlist(lapply(components,function(x){x==")"})))
components <- components[[1]][(locstart+1):(locend-1)]
mn <- length(components)
if (mn == 4){
components <- c(components[1],paste0(components[2],components[3]),components[mn])
}
ttl <- "ES components"
# Check for xreg
tempInit <- x[[AL,which(colnames(x)=="initial")]]
x.n[AL] <- length(tempInit$xreg)
}
# Error term
if (components[1]=="A"){
comps[AL,1] <- 1
} else {
comps[AL,1] <- 2
}
comps.char[AL,1] <- components[1]
# Trend term
if (components[2]=="A"){
comps[AL,2] <- 1
} else if (components[2]=="Ad"){
comps[AL,2] <- 1.5
} else if (components[2]=="M"){
comps[AL,2] <- 2
} else if (components[2]=="Md"){
comps[AL,2] <- 2.5
} else {
comps[AL,2] <- 0
}
comps.char[AL,2] <- components[2]
# Season term
if (components[3]=="A"){
comps[AL,3] <- 1
} else {if (components[3]=="M"){
comps[AL,3] <- 2
} else
comps[AL,3] <- 0
}
comps.char[AL,3] <- components[3]
comps[AL,4] <- x[[AL,idx.AL]]
}
# Add row for xreg if needed
if (sum(x.n)>0 & xreg.plot == TRUE){
x.add <- 1
comps <- cbind(comps[,1:3],x.n>0,comps[,4])
} else {
x.add <- 0
}
image(min(comps[,4+x.add]):max(comps[,4+x.add]), 1:(3+x.add), matrix(comps[,(3+x.add):1],ncol=(3+x.add)), axes=FALSE, col=brewer.pal(8,"PuBu")[2:6], #rev(heat.colors(5)),
ylab="Components", xlab="Aggregation Level", main=ttl,breaks=c(-1,0,1,1.5,2,2.5))
axis(2, at=1:(3+x.add), labels=list("Error","Trend","Season","Xreg")[(3+x.add):1])
axis(1, at=min(comps[,4+x.add]):max(comps[,4+x.add]))
# Grey seasonality
k <- 1+x.add
for (AL in which(!perm_seas)){
polygon(c((AL-0.5)*c(1,1),(AL+0.5)*c(1,1)),c(.5+k,-.5+k,-.5+k,.5+k),col="gray60",border=NA)
}
box()
for (i in 1:(4+x.add)){
for (AL in 1:max(ALplot)){
if (i==1){
lines(c(AL-0.5+minimumAL-1,AL-0.5+minimumAL-1),c(0,4+x.add),col="black")
}
if (i<4 & AL<=max(comps[,4+x.add])){
text(AL+minimumAL-1,i+x.add,comps.char[AL,4-i])
}
if (x.add==1){
text(AL+minimumAL-1,1,x.n[AL])
}
}
lines(c(min(comps[,4+x.add])-0.5,max(comps[,4+x.add])+0.5),c(i-0.5,i-0.5),col="black")
}
}
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.