Piecewise linear trends

October 27, 2015
By

[This article was first published on Hyndsight » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



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)