Site icon R-bloggers

Piecewise linear trends

[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:

   

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

If the associated coefficients of and are and , then gives the slope of the trend before time , while the slope of the line after time is given by .

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

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 .

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 , we can set

   

xreg <- pmin(0, x1-tau)

where the first 0 in the column is in row . Additional knots can be handled in the same way. For example, if there are two knots, then will be the slope of the trend up to the first knot, and 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.