Estimating a nonlinear time series model in R


21 January 2014

time series

There are quite a few R packages available for nonlinear time series analysis, but sometimes you need to code your own models. Here is a simple example to show how it can be done.

The model is a first order threshold autoregression:

y_t = \begin{cases} \alpha y_{t-1} + e_t & \text{if $y_{t-1} \le r$}\\ \beta y_{t-1} + \gamma e_t & \text{if $y_{t-1}>r$} \end{cases}

where e_t is a Gaussian white noise series with variance \sigma^2. The following function will generate some random data from this model.

simnlts <- function(n, alpha, beta, r, sigma, gamma, burnin=100)
  # Generate noise
  e <- rnorm(n+burnin, 0, sigma)
  # Create space for y
  y <- numeric(n+burnin)
  # Generate time series
  for(i in 2:(n+burnin))
    if(y[i-1] <= r)
      y[i] <- alpha*y[i-1] + e[i]
      y[i] <- beta*y[i-1] + gamma*e[i]
  # Throw away first burnin values
  y <- ts(y[-(1:burnin)])
  # Return result

The “burn-in” parameter allows the first value of the series to be a random draw from the stationary distribution of the process (assuming it is stationary).

We will estimate the model by minimizing the sum of squared errors across both regimes. Other approaches that take account of the different variances in the two regimes may be better, but it is useful to keep it simple, at least initially. The following function uses a non-linear optimizer to estimate \alpha, \beta and r. To ensure good starting values, we begin the optimization with \alpha=\beta=0 and set the initial value of r to be the mean of the observed data.

fitnlts <- function(x)
  ss <- function(par,x)
    alpha <- par[1]
    beta <- par[2]
    r <- par[3]
    n <- length(x)
    # Check that each regime has at least 10% of observations
    if(sum(x<=r) < n/10 | sum(x>r) < n/10)
    e1 <- x[2:n] - alpha*x[1:(n-1)]
    e2 <- x[2:n] - beta*x[1:(n-1)]
    regime1 <- (x[1:(n-1)] <= r)
    e <- e1*(regime1) + e2*(!regime1)
  fit <- optim(c(0,0,mean(x)),ss,x=x,control=list(maxit=1000))
  if(fit$convergence > 0)

There are a few things to notice here.

The functions can be used as follows (with some example parameters):

y <- simnlts(100, 0.5, -1.8, -1, 1, 2)

A similar approach can be used for other time series models.