# R/idclass.R In tsintermittent: Intermittent Time Series Forecasting

#### Documented in idclass

```idclass <- function(data,type=c("PKa","SBC","KHa","KH","PK"),a.in=NULL,
outplot=c("summary","detail","none"),plot.focus=NULL){
# Time series categorisation for intermittent demand
#
# Inputs
#   data          Time series dataset. Each column is a series.
#                 Alternatively this can be a single series.
#   type          Type of categorisation:
#                   "SBC" - Syntetos Boylan Croston;
#                   "KH"  - Kostenko Hyndman (exact*);
#                   "KHa" - Kostenko Hyndman (approximate);
#                   "PK"  - Petropoulos Kourentzes (exact*);
#                   "PKa" - Petropoulos Kourentzes (approximate).
#                 * These are computationally expensive, as SBA is optimised
#                   for each time series.
#   a.in          Vector of SBA demand interval smoothing parameters. This
#                 must be same length as number of series. This is used for
#                 categorisations "KH" and "PK". If a.in == NULL then the
#                 parameters are calculated internally using MAR as a cost
#                 function.
#   outplot       Plot results of categorisation:
#                   "summary" - simlified plot that reports number of series
#                               in each class and cut-off points;
#                   "detail"  - scatterplot between average inter-demand
#                               interval (p) and squared coefficient of
#                               variation of non-zero demand (CV^2). Series
#                               that are categorised for SBA or SES are plotted
#                   "none"    - do not produce plot.
#   plot.focus    Only relevant to outplot == "detail". Can be used to specify
#                 the maximum p and CV^2 to plot, so that the scatterplot can
#                 be focused on the separation area between the categories. Use
#                 vector of two elements. First one is max p and second one is
#                 max CV^2. Example: plot.focus=c(1.5,1.5). If NULL then maximums
#                 are defined from the dataset.
#
# Outputs
#   idx.croston   Index of series that are categorised under Croston.
#   idx.sba       Index of series that are categorised under SBA.
#   idx.ses       Index of series that are categorised under SES. Provided only
#                 for "PK" and "PKa" types.
#   cv2           Coefficient of variation squared of non-zero demand.
#   p             Inter-demand interval.
#   summary       Summary of number of series under each category.
#
# Example:
#   # Create/load some data. Each column is a time series
#   dataset <- simID(100,60,idi=1.15,cv2=0.3)
#   idclass(dataset)
#
# Notes:
# Classification schemes described in:
# F. Petropoulos and N. Kourentzes, 2014, Journal of Operational Research Society.
# http://dx.doi.org/10.1057/jors.2014.62
# http://kourentzes.com/forecasting/2014/05/13/forecast-combinations-for-intermittent-demand/
#
# Nikolaos Kourentzes, 2014 <nikolaos@kourentzes.com>

type <- type
type <- toupper(type)
outplot <- outplot

# Initialise
if ((sum(dim(data)==1)>0 | is.vector(data)) && !isa(data,"data.frame")){
data <- matrix(data,nrow=length(data))
N <- 1
p <- NA
v <- NA
if (is.null(a.in)){
a <- NA
} else {
a <- a.in
}
} else {
N <- dim(data)
p <- vector("numeric",N)    # Average interval
v <- vector("numeric",N)    # Squared CV
if (is.null(a.in)){
a <- vector("numeric",N)  # Interval alpha with SBA
} else {
a <- a.in
if (length(a)!=N){
stop("Length of a.in must be equal to number of series.")
}
}
}

for (i in 1:N){
# Get time series
temp <- data[,i]
temp <- temp[!is.na(temp)]
# Croston decomposition
nzd <- which(temp!=0)
k <- length(nzd)
z <- temp[nzd]                        # Demand
x <- c(nzd,nzd[2:k]-nzd[1:(k-1)])  # Intervals
p[i] <- mean(x)
v[i] <- (sd(z)/mean(z))^2
if (is.null(a.in)){
if (type=="KH" | type=="PK"){
fit <- crost(temp,h=0,type="sba")
a[i] <- fit\$weights
}
}
}

# Classify between SBA & Croston
if (type=="KH" | type=="PK"){
# Exact Konstenko & Hyndman
use.sba <- v > (4*p*(2-p)-a*(4-a)-p*(p-1)*(4-a)*(2-a))/(p*(4-a)*(2*p-a))
} else if(type=="SBC"){
# Syntetos, Boylan & Croston
use.sba <- (p > 1.32 | v > 0.49)
} else {
# Approximate Konstenko & Hyndman
use.sba <- v > (2-(3/2)*p)
}
use.croston <- !use.sba

# Classify between SES and Croston/SBA
if (type=="PK" | type=="PKA"){
use.ses <- p <= 1
use.sba[use.ses] <- FALSE
use.croston[use.ses] <- FALSE
idx.ses <- which(use.ses)
} else {
use.ses <- NULL
idx.ses <- NULL
}
idx.croston <- which(use.croston)
idx.sba <- which(use.sba)

# Summarise results
if (is.null(use.ses)){
summary <- rbind(sum(use.croston),sum(use.sba))
rownames(summary) <- c("Croston","SBA")
} else {
summary <- rbind(sum(use.croston),sum(use.sba),sum(use.ses))
rownames(summary) <- c("Croston","SBA","SES")
}
colnames(summary) <- "Series"

if (type=="PKA"){
type.plot = "PKa"
} else if(type=="KHA"){
type.plot = "KHa"
} else {
type.plot = type
}

# Produce detailed plots
if (outplot == "detail"){
# Set space in plot for PK and PKa types
if (type=="PK" | type=="PKA"){
if (!is.null(plot.focus)){
xmin <- 1-(plot.focus-1)*0.1
} else {
xmin <- 1-(max(c(max(p)+diff(range(p))*0.1),1.5)-1)*0.1
}
} else {
xmin <- 1
}
# Set limits depending if plot.focus is given or not
if (!is.null(plot.focus)){
xx <- c(xmin,plot.focus)
yy <- c(0,plot.focus)
plot.out <- (p > plot.focus | v > plot.focus)
} else {
xx <- c(xmin,max(c(max(p)+diff(range(p))*0.1),1.5))
yy <- c(0,max(c(max(v)+diff(range(v))*0.1),1.5))
plot.out <- vector("logical",N)
}
# Plot SBA series in focus
plot(p[use.sba & !plot.out],v[use.sba & !plot.out],
type="p",pch=20,xaxs="i",yaxs="i",
xlim=xx,ylim=yy,cex=0.8,xlab="p",ylab=parse(text="CV^2"),
main=paste(type.plot,"classification"))
# Plot SBA series out-of-focus
if (sum(plot.out)>1){
p.temp <- p[plot.out]
p.temp[p.temp>xx] <- xx
v.temp <- v[plot.out]
v.temp[v.temp>yy] <- yy
points(p.temp,v.temp,pch=21,bg="darkgray",cex=1)
}
# Create shaded area for Croston series
if (type == "SBC"){
polygon(c(1,1.32,1.32,1),c(0,0,0.49,0.49),col="lightgrey",border=NA)
} else if(type == "KHA" | type == "PKA"){
polygon(c(1,4/3,1),c(0,0,0.5),col="lightgrey",border=NA)
} else if(type == "KH" | type == "PK"){
p.x <- c(seq(1,4/3,0.02),4/3)
a.p <- 1
v.y <- (4*p.x*(2-p.x)-a.p*(4-a.p)-p.x*(p.x-1)*(4-a.p)*(2-a.p))/(p.x*(4-a.p)*(2*p.x-a.p))
polygon(c(p.x,1),c(v.y,0),col="lightgrey",border=NA)
a.p <- 0
v.y2 <- (4*p.x*(2-p.x)-a.p*(4-a.p)-p.x*(p.x-1)*(4-a.p)*(2-a.p))/(p.x*(4-a.p)*(2*p.x-a.p))
v.y2[v.y2<0] <- 0
polygon(c(p.x,p.x[18:1]),c(v.y,v.y2[18:1]),col=gray(0.2,0.3),border=NA)
}
# Plot Croston series
points(p[use.croston],v[use.croston],pch=1,col="red",cex=0.6)
# Create area for ses
if (type == "PKA" | type == "PK"){
polygon(c(0,1,1,0),c(0,0,yy,yy),col=gray(0.4),border=NA)
points(rep(xmin+(1-xmin)/2,sum(use.ses & !plot.out)),v[use.ses & !plot.out],
col="white",pch=20,cex=0.8)
if (sum(use.ses & plot.out)>1){
points(xmin+(1-xmin)/2,yy,col="grey",pch=21,bg="white")
}
}
}
# Produce simplified plots
if (outplot=="summary"){
# Set plot area
if (type=="PK" | type=="PKA"){
xx <- c(0.94,1.5)
} else if(type=="SBC"){
xx <- c(1,1.64)
} else {
xx <- c(1,1.5)
}
yy <- c(0,1)
# Start a plot
plot(0,0,xlim=xx,ylim=yy,xaxs="i",yaxs="i",xlab="p",ylab=parse(text="CV^2"),
main=paste(type.plot,"classification"),xaxt="n",yaxt="n")
# Create shaded area for Croston series
if (type == "SBC"){
polygon(c(1,1.32,1.32,1),c(0,0,0.49,0.49),col="lightgrey",border=NA)
} else if(type == "KHA" | type == "PKA"){
polygon(c(1,4/3,1),c(0,0,0.5),col="lightgrey",border=NA)
} else if(type == "KH" | type == "PK"){
p.x <- c(seq(1,4/3,0.02),4/3)
a.p <- 1
v.y <- (4*p.x*(2-p.x)-a.p*(4-a.p)-p.x*(p.x-1)*(4-a.p)*(2-a.p))/(p.x*(4-a.p)*(2*p.x-a.p))
polygon(c(p.x,1),c(v.y,0),col="lightgrey",border=NA)
a.p <- 0
v.y2 <- (4*p.x*(2-p.x)-a.p*(4-a.p)-p.x*(p.x-1)*(4-a.p)*(2-a.p))/(p.x*(4-a.p)*(2*p.x-a.p))
v.y2[v.y2<0] <- 0
polygon(c(p.x,p.x[18:1]),c(v.y,v.y2[18:1]),col=gray(0.2,0.3),border=NA)
}
# Label ticks
if (type != "SBC"){
axis(1,at=c(1,4/3),labels=c("1","4/3"))
axis(2,at=c(0,0.5),labels=c("0","1/2"))
} else {
# Special plot for SBC
axis(1,at=c(1,1.32),labels=c("1","1.32"))
axis(2,at=c(0,0.49),labels=c("0","0.49"))
lines(c(1,xx),c(0.49,0.49),lty=3,col="grey")
lines(c(1.32,1.32),c(0,yy),lty=3,col="grey")
text(1.16,0.75,"Erratic",col=gray(0.4))
text(1.16,0.25,"Smooth",col=gray(0.4))
text(1.48,0.75,"Lumpy",col=gray(0.4))
text(1.48,0.25,"Intermittent",col=gray(0.4))
}
# Report number of series in each area
# Shaded area and text for PK and PKa classifications
if (type == "PKA" | type == "PK"){
polygon(c(0,1,1,0),c(0,0,yy,yy),col=gray(0.6),border=NA)
lines(c(0.94,1),c(0.5,0.5),lty=3)
}
}

return(list(idx.croston=idx.croston,idx.sba=idx.sba,idx.ses=idx.ses,
cv2=v,p=p,summary=summary))

}
```

## Try the tsintermittent package in your browser

Any scripts or data that you put into this service are public.

tsintermittent documentation built on July 18, 2022, 9:06 a.m.