**R on Rob J Hyndman**, 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

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

where \(\boldsymbol{y}=[y_1,\dots,y_T]'\), \(\boldsymbol{n}=[n_1,\dots,n_T]'\), \(\boldsymbol{\beta}=[\beta_1]\) and \(\boldsymbol{X} = [x_{1,1},\dots,x_{1,T}]'\). Note that I have left the intercept \(\beta\_0\) out of the vector \(\boldsymbol{\beta}\) so that the \(\boldsymbol{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:

\[

\boldsymbol{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 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.

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

\[

\boldsymbol{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')
```

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 \(\boldsymbol{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

\[

\boldsymbol{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.

**leave a comment**for the author, please follow the link and comment on their blog:

**R on Rob J Hyndman**.

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...