Piecewise linear trends

October 27, 2015
By

(This article was first published on Hyndsight » R, and kindly contributed to R-bloggers)

I prepared the following notes for a consulting client, and I thought they might be of interest to some other people too.

Let y_t denote the value of the time series at time t, and suppose we wish to fit a trend with correlated errors of the form

    \[y_t = f(t) + n_t,\]

where f(t) represents the possibly nonlinear trend and n_t is an autocorrelated error process.

For example, if f(t) = \beta_0+\beta_1 t is a linear function, then we can simply set x_{1,t}=t and define

    \[y_t = \beta_0 + \beta_1x_{1,t} + n_t.\]

In matrix form we can write

    \[\bm{y} = \beta_0 + \bm{X}\bm{\beta} + \bm{n},\]

where \bm{y}=[y_1,\dots,y_T]', \bm{n}=[n_1,\dots,n_T]', \bm{\beta}=[\beta_1] and \bm{X} = [x_{1,1},\dots,x_{1,T}]'. Note that I have left the intercept \beta_0 out of the vector \bm{\beta} so that the \bm{X} matrix matches the required xreg argument in auto.arima.

This model can be estimated by setting the xreg argument to be a matrix with one column:

    \[\bm{X} = \left[\begin{array}{c} 1\\ 2\\ 3\\ 4\\ \vdots\\ T \end{array}\right]\]

x1 <- 1:length(y)
fit <- auto.arima(y, xreg=x1)

The associated coefficient is the slope of the trend line.

Here is a simple example of a linear trend fitted to the Asian sheep data from the fpp package :

library(fpp)
T <- length(livestock)
x1 <- seq(T)
fit <- auto.arima(livestock, xreg=x1)
fc <- forecast(fit, xreg=T+seq(10))
b0 <- coef(fit)["intercept"]
b1 <- coef(fit)["x1"]
t <- seq(T+10)
trend <- ts(b0 + b1*t, start=start(livestock))
 
plot(fc, main="Linear trend with AR(1) errors")
lines(trend, col='red')
A linear trend fitted to the Asian sheep data. The automatically selected error term is an AR(1) process.

A linear trend fitted to the Asian sheep data. The automatically selected error term is an AR(1) process.

A more flexible approach is to use a piecewise linear trend which bends at some time. If the trend bends at time \tau, then it can be specified by including the following predictors in the model.

(1)   \begin{align*}  x_{1,t}  &=  t \\  x_{2,t} &=  \begin{cases} 0 & t < \tau;\\ (t-\tau) &  t \ge \tau. \end{cases} \end{align*}

In auto.arima, set xreg to be a matrix with two columns:

    \[\bm{X} = \left[\begin{array}{ll} 1 & 0\\ 2 & 0\\ 3 & 0\\ 4 & 0\\ \vdots\\ \tau & 0 \\ \tau+1 & 1\\ \tau+2 & 2\\ \vdots \\ T & T-\tau \end{array}\right]\]

fit <- auto.arima(y, xreg=cbind(x1, pmax(0,x1-tau))

If the associated coefficients of x_{1,t} and x_{2,t} are \beta_1 and \beta_2, then \beta_1 gives the slope of the trend before time \tau, while the slope of the line after time \tau is given by \beta_1+\beta_2.

This can be extended to allow any number of “bend points” known as knots. Just add additional columns with 0s before each knot, and values 1, 2, … after the knot.

Here is a piecewise linear trend fitted to the Asian sheep data with knots at years 1990 and 1992:

x2 <- pmax(0, x1-30)
x3 <- pmax(0, x1-32)
fit <- auto.arima(livestock, xreg=cbind(x1,x2,x3))
fc <- forecast(fit, 
   xreg=cbind(max(x1)+seq(10), max(x2)+seq(10), max(x3)+seq(10)))
b0 <- coef(fit)["intercept"]
b1 <- coef(fit)["x1"]
b2 <- coef(fit)["x2"]
b3 <- coef(fit)["x3"]
trend <- ts(b0 + b1*t + b2*pmax(0,t-30) + b3*pmax(0,t-32),
    start=start(livestock))
 
plot(fc, main="Piecewise linear trend with AR(1) errors")
lines(trend, col='red')
A piecewise-linear trend fitted to the Asian sheep data

A piecewise-linear trend fitted to the Asian sheep data

If there is to be no trend before the first knot, but a piecewise linear trend thereafter, leave out the first column of the above matrix \bm{X}.

If there is to be a piecewise linear trend up to the last knot, but no trend thereafter, a slightly modified set up can be used. For one knot at time \tau, we can set

    \[\bm{X} = \left[\begin{array}{r} 1-\tau \\ 2-\tau \\ \vdots\\ -2\\ -1\\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right]\]

xreg <- pmin(0, x1-tau)

where the first 0 in the column is in row \tau. Additional knots can be handled in the same way. For example, if there are two knots, then \beta_1+\beta_2 will be the slope of the trend up to the first knot, and \beta_2 will be the slope between the first and second knots.

To leave a comment for the author, please follow the link and comment on their blog: Hyndsight » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)