My new book is Forecasting: principles and practice, co-authored with George Athanasopoulos. It is available online and free-of-charge. We have written about 2/3 of the book so far (all of which is already available online), and we plan to finish it by the end of 2012. We hope to make a print version of the book available on Amazon in early 2013.

This textbook is intended to provide a comprehensive introduction to forecasting methods and present enough information about each method for readers to use them sensibly. We don’t attempt to give a thorough discussion of the theoretical details behind each method, although the references at the end of each chapter will fill in many of those details. We use R throughout the book and we intend students to learn how to forecast with R.

The book has it’s own R package: fpp. This contains all the data sets used in the book, and also loads a few other packages that are necessary to complete the examples.

The book is different from other forecasting textbooks in several ways.

• It is free and online, making it accessible to a wide audience.
• It is based around the forecast package for R.
• It is continuously updated. You don’t have to wait until the next edition for errors to be removed or new methods to be discussed. We will update the book frequently.
• There are dozens of real data examples taken from our own consulting practice. We have worked with hundreds of businesses and organizations helping them with forecasting issues, and this experience has contributed directly to many of the examples given here, as well as guiding our general philosophy of forecasting.
• We emphasise graphical methods more than most forecasters. We use graphs to explore the data, analyse the validity of the models fitted and present the forecasting results.

There are now several blog aggregators available that might be of interest to readers here. And this blog is now syndicated on several other sites including those listed below.

• R-bloggers: for all R-related blogs. The posts tagged R from this blog are syndicated there along with about 300 other R blogs.
• Statsblogs for statistical blogs. This is a very new aggregator, but is growing fast. There is naturally some overlap with R-bloggers. All posts from this blog are syndicated there.
• TeX community for TeX related blogs. The posts tagged LaTeX from this blog are syndicated there along with about 40 other TeX blogs.
• Mathblogging.org aggregates a number of mathematics blogs. All posts from this blog are syndicated there.

In addition, for those interested in economics, EconAcademics.org aggregates a number of economics blogs (but does not include this blog as I rarely post anything about economics).

I have also set up two aggregations for new research papers:

• StatisticsPapers includes all papers appearing in over 60 of the major statistics journals as well as the statistics section of arXiv.
• Forecasting papers includes all new papers on forecasting that have appeared on RePEc or in either the International Journal of Forecasting or the Journal of Forecasting. The papers from RePEc constitute the weekly NEP-FOR report.

So if all the blogs around are overwhelming, and you don’t know where to start, select one or two of these aggregators and you’re off and running.

If you’ve no idea how to subscribe to a blog, see my post on using Google Reader.

If you want help, please ask at either stats.stackexchange.com (for R or statistics questions) or tex.stackexchange.com (for LaTeX questions).

Unless you are one of my students, the only questions I will answer are ones that concern my R packages or research papers. And even then, I won’t reply if the answer is in the help files. I write those help files for a reason, so please read them.

I’m sorry I can’t do more, but if I did everything people ask me to do, I’d never write any papers or produce any R packages, and I think that’s a better use of my time.

1. Wang, Smith and Hyndman (2006) Characteristic-​​based clustering for time series data. Data Mining and Knowledge Discovery, 13(3), 335–364.
2. Wang, Smith-Miles and Hyndman (2009) “Rule induction for forecasting method selection: meta-​​learning the characteristics of univariate time series”, Neurocomputing, 72, 2581–2594.

I’ve since had a lot of requests for the code which one of my coauthors has been helpfully emailing to anyone who asked. But to make it easier, we thought it might be helpful if I post some updated code here. This is not the same as the R code we used in the paper, as I’ve improved it in several ways (so it will give different results). If you just want the code, skip to the bottom of the post.

Finding the period of the data

Usually in time series work, we know the period of the data (if the observations are monthly, the period is 12, for example). But in this project, some of our data was of unknown period and we wanted a method to automatically determine the appropriate period. The method we used was based on local peaks and troughs in the ACF. But I’ve since devised a better approach (prompted on crossvalidated.com) using an estimate of the spectral density:

find.freq <- function(x)
{
  n <- length(x)
  spec <- spec.ar(c(na.contiguous(x)),plot=FALSE)
  if(max(spec$spec)>10) # Arbitrary threshold chosen by trial and error.
  {
    period <- round(1/spec$freq[which.max(spec$spec)])
    if(period==Inf) # Find next local maximum
    {
      j <- which(diff(spec$spec)>0)
      if(length(j)>0)
      {
        nextmax <- j[1] + which.max(spec$spec[j[1]:500])
        if(nextmax <= length(spec$freq))
          period <- round(1/spec$freq[nextmax])
        else
          period <- 1
      }
      else
        period <- 1
    }
  }
  else
    period <- 1
 
  return(period)
}

 
The function is called find.freq because time series people often call the period of seasonality the “frequency” (which is of course highly confusing).

Decomposing the data into trend and seasonal components

We needed a measure of the strength of trend and the strength of seasonality, and to do this we decomposed the data into trend, seasonal and error terms.

Because not all data could be decomposed additively, we first needed to apply an automated Box-Cox transformation. We tried a range of Box-Cox parameters on a grid, and selected the one which gave the most normal errors. That worked ok, but I’ve since found some papers that provide quite good automated Box-Cox algorithms that I’ve implemented in the forecast package. So this code uses Guerrero’s (1993) method instead.

For seasonal time series, we decomposed the transformed data using an stl decomposition with periodic seasonality.

For non-seasonal time series, we estimated the trend of the transformed data using penalized regression splines via the mgcv package.

decomp <- function(x,transform=TRUE)
{
  require(forecast)
  # Transform series
  if(transform & min(x,na.rm=TRUE) >= 0)
  {
    lambda <- BoxCox.lambda(na.contiguous(x))
    x <- BoxCox(x,lambda)
  }
  else
  {
    lambda <- NULL
    transform <- FALSE
  }
  # Seasonal data
  if(frequency(x)>1)
  {
    x.stl <- stl(x,s.window="periodic",na.action=na.contiguous)
    trend <- x.stl$time.series[,2]
    season <- x.stl$time.series[,1]
    remainder <- x - trend - season
  }
  else #Nonseasonal data
  {
    require(mgcv)
    tt <- 1:length(x)
    trend <- rep(NA,length(x))
    trend[!is.na(x)] <- fitted(gam(x ~ s(tt)))
    season <- NULL
    remainder <- x - trend
  }
  return(list(x=x,trend=trend,season=season,remainder=remainder,
    transform=transform,lambda=lambda))
}

Putting everything on a [0,1] scale

We wanted to measure a range of characteristics such as strength of seasonality, strength of trend, level of nonlinearity, skewness, kurtosis, serial correlatedness, self-similarity, level of chaoticity (is that a word?) and the periodicity of the data. But we wanted all these on the same scale which meant mapping the natural range of each measure onto [0,1]. The following two functions were used to do this.

# f1 maps [0,infinity) to [0,1]
f1 <- function(x,a,b)
{
  eax <- exp(a*x)
  if (eax == Inf)
    f1eax <- 1
  else
    f1eax <- (eax-1)/(eax+b)
  return(f1eax)
}
 
# f2 maps [0,1] onto [0,1]
f2 <- function(x,a,b)
{
  eax <- exp(a*x)
  ea <- exp(a)
  return((eax-1)/(eax+b)*(ea+b)/(ea-1))
}

 
The values of and in each function were chosen so the measure had a 90th percentile of 0.10 when the data were iid standard normal, and a value of 0.9 using a well-known benchmark time series.

Calculating the measures

Now we are ready to calculate the measures on the original data, as well as on the adjusted data (after removing trend and seasonality).

measures <- function(x)
{
  require(forecast)
 
  N <- length(x)
  freq <- find.freq(x)
  fx <- c(frequency=(exp((freq-1)/50)-1)/(1+exp((freq-1)/50)))
  x <- ts(x,f=freq)
 
  # Decomposition
  decomp.x <- decomp(x)
 
  # Adjust data
  if(freq > 1)
    fits <- decomp.x$trend + decomp.x$season
  else # Nonseasonal data
    fits <- decomp.x$trend
  adj.x <- decomp.x$x - fits + mean(decomp.x$trend, na.rm=TRUE)
 
  # Backtransformation of adjusted data
  if(decomp.x$transform)
    tadj.x <- InvBoxCox(adj.x,decomp.x$lambda)
  else
    tadj.x <- adj.x
 
  # Trend and seasonal measures
  v.adj <- var(adj.x, na.rm=TRUE)
  if(freq > 1)
  {
    detrend <- decomp.x$x - decomp.x$trend
    deseason <- decomp.x$x - decomp.x$season
    trend <- ifelse(var(deseason,na.rm=TRUE) < 1e-10, 0, 
      max(0,min(1,1-v.adj/var(deseason,na.rm=TRUE))))
    season <- ifelse(var(detrend,na.rm=TRUE) < 1e-10, 0,
      max(0,min(1,1-v.adj/var(detrend,na.rm=TRUE))))
  }
  else #Nonseasonal data
  {
    trend <- ifelse(var(decomp.x$x,na.rm=TRUE) < 1e-10, 0,
      max(0,min(1,1-v.adj/var(decomp.x$x,na.rm=TRUE))))
    season <- 0
  }
 
  m <- c(fx,trend,season)
 
  # Measures on original data
  xbar <- mean(x,na.rm=TRUE)
  s <- sd(x,na.rm=TRUE)
 
  # Serial correlation
  Q <- Box.test(x,lag=10)$statistic/(N*10)
  fQ <- f2(Q,7.53,0.103)
 
  # Nonlinearity
  p <- terasvirta.test(na.contiguous(x))$statistic
  fp <- f1(p,0.069,2.304)
 
  # Skewness
  sk <- abs(mean((x-xbar)^3,na.rm=TRUE)/s^3)
  fs <- f1(sk,1.510,5.993)
 
  # Kurtosis
  k <- mean((x-xbar)^4,na.rm=TRUE)/s^4
  fk <- f1(k,2.273,11567)
 
  # Hurst=d+0.5 where d is fractional difference.
  H <- fracdiff(na.contiguous(x),0,0)$d + 0.5
 
  # Lyapunov Exponent
  if(freq > N-10)
      stop("Insufficient data")
  Ly <- numeric(N-freq)
  for(i in 1:(N-freq))
  {
    idx <- order(abs(x[i] - x))
    idx <- idx[idx < (N-freq)]
    j <- idx[2]
    Ly[i] <- log(abs((x[i+freq] - x[j+freq])/(x[i]-x[j])))/freq
    if(is.na(Ly[i]) | Ly[i]==Inf | Ly[i]==-Inf)
      Ly[i] <- NA
  }
  Lyap <- mean(Ly,na.rm=TRUE)
  fLyap <- exp(Lyap)/(1+exp(Lyap))
 
  m <- c(m,fQ,fp,fs,fk,H,fLyap)
 
  # Measures on adjusted data
  xbar <- mean(tadj.x, na.rm=TRUE)
  s <- sd(tadj.x, na.rm=TRUE)
 
  # Serial
  Q <- Box.test(adj.x,lag=10)$statistic/(N*10)
  fQ <- f2(Q,7.53,0.103)
 
  # Nonlinearity
  p <- terasvirta.test(na.contiguous(adj.x))$statistic
  fp <- f1(p,0.069,2.304)
 
  # Skewness
  sk <- abs(mean((tadj.x-xbar)^3,na.rm=TRUE)/s^3)
  fs <- f1(sk,1.510,5.993)
 
  # Kurtosis
  k <- mean((tadj.x-xbar)^4,na.rm=TRUE)/s^4
  fk <- f1(k,2.273,11567)
 
  m <- c(m,fQ,fp,fs,fk)
  names(m) <- c("frequency", "trend","seasonal",
    "autocorrelation","non-linear","skewness","kurtosis",
    "Hurst","Lyapunov",
    "dc autocorrelation","dc non-linear","dc skewness","dc kurtosis")
 
  return(m)
}

Here is a quick example applied to Australian monthly gas production:

library(forecast)
measures(gas)
     frequency              trend           seasonal    autocorrelation 
        0.1096             0.9989             0.9337             0.9985 
    non-linear           skewness           kurtosis              Hurst 
        0.4947             0.1282             0.0055             0.9996 
      Lyapunov dc autocorrelation      dc non-linear        dc skewness 
        0.5662             0.1140             0.0538             0.1743 
   dc kurtosis 
        0.9992

 
The function is far from perfect, and it is not hard to find examples where it fails. For example, it doesn’t work with multiple seasonality — try measure(taylor) and check the seasonality. Also, I’m not convinced the kurtosis provides anything useful here, or that the skewness measure is done in the best way possible. But it was really a proof of concept, so we will leave it to others to revise and improve the code.

In our papers, we took the measures obtained using R, and produced self-organizing maps using Viscovery. There is now a som package in R for that, so it might be possible to integrate that step into R as well. The dendogram was generated in matlab, although that could now also be done in R using the ggdendro package for example.

Download the code in a single file.

Now there is a great new website for LaTeX templates: www.latextemplates.com. There are some nice templates for letters, lab reports, calendars, theses, assignments, essays, and CVs.  The templates are well-structured with lots of comments to make it easy to understand how they work, and to make modifications. Even experienced LaTeXers will probably learn some new tricks and new packages from browsing the templates.

For example, a search on “forecasting” yields the following results:

The h-index is the largest number such that at least articles in that publication were cited at least times each. The h5-index is the h-index calculated using only articles published in the last five complete calendar years (2007–2011 in this example). So there were 29 articles published in the IJF between 2007 and 2011 that have each been cited at least 29 times.

The h5-median is the median number of citations of the articles making up the h5-index. That is, for the IJF the 29 articles have received a median 43 citations. By clicking on the number in the h5-index column, you can see which articles have been included.

Currently, any discussion of journal metrics is dominated by the ISI 2-year impact factor, equal to the average number of citations received per paper published in that journal during the two preceding years. In my view the h5-index is a far better measure than the 2-year IF. Here are some reasons why.

1. I like the fact that a 5-year index has been used, rather than the 2-year impact factor favoured by the ISI. Two years is much too short, and leads to a lot of year-to-year variation, at least in the fields I’m interested in. Five years provides a smoother measure that will not change so much from year to year, yet will still represent recent quality rather than what the journal might have been like many years ago.
2. The h5-index cannot be dominated by one star paper. If there is only one great paper in the journal, and everything else is not cited at all, the h5-index will be 1. On the other hand, the ISI 2-year impact factor will be greatly increased by that single paper due to the non-robustness of the mean.
3. It is harder to game the h5-index than the ISI impact factor. One of the games that some journals play is to force authors wanting their article published in the journal to cite other articles published in the journal, thus artificially increasing the number of citations. This has a direct and immediate impact on the ISI impact factor as it is based on average citations per article. But it will be harder to use this game on the journal h5-index.To see this, imagine that the Journal of Forecasting wanted to improve their h5-index to 30 and so beat the IJF. They currently only have 5 papers with 30 or more citations, so they would need to get another 25 papers up to that level, 14 of which currently have fewer than 16 citations. So that’s more than 14 extra JF citations for each of those 14 papers — even with bogus citations that’s not going to happen. (I am not suggesting that JF ever plays such games, only pointing out that it will be much harder for journals to game the h5-index than the ISI 2-year impact factor.)

In summary, the h5-index is simple to understand, hard to manipulate, and provides a reasonable if crude measure of the respect accorded to a journal by scholars within its field. While journal metrics are no guarantee of the quality of a journal, if they are going to be used we should use the best available, and Google’s h5-index is a big improvement on the ISI impact factor.

In the following two posts, Frank Davenport shows how it can be done:

1. Plotting forecast() objects in ggplot part 1: Extracting the Data
2. Plotting forecast() objects in ggplot part 2: Visualize Observations, Fits, and Forecasts

Fortunately, it is not difficult. Go to Options/Commands where all the commands used by TeXstudio are specified. You probably don’t need standard LaTeX these days, so replace the LaTeX command with the following.

xelatex -interaction=nonstopmode %.tex

Then click OK.

Now the LaTeX button at the top of the screen is mapped to XeLaTeX rather than standard LaTeX. You can still access pdfLaTeX via its button for your non-XeLaTeX files.

Before anyone comments that you need standard LaTeX for when eps graphics are used, see Converting eps to pdf.

The classics are also worth reading, and remain relevant despite the 20 or 30 years that have elapsed since they appeared.