Статьи

Оценка модели нелинейного временного ряда в R

Март 31, 2018

Существует немало пакетов R, доступных для нелинейного анализа временных рядов, но иногда вам нужно кодировать собственные модели. Вот простой пример, чтобы показать, как это можно сделать.

Модель представляет собой пороговую авторегрессию первого порядка:

где e_tгауссовский ряд белого шума с дисперсией \ Sigma ^ 2. Следующая функция будет генерировать случайные данные из этой модели.

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]
    else
      y[i] <- beta*y[i-1] + gamma*e[i]
  }
  # Throw away first burnin values
  y <- ts(y[-(1:burnin)])
  # Return result
  return(y)
}

Параметр «прожигания» позволяет первому значению ряда быть случайным извлечением из стационарного распределения процесса (при условии, что оно стационарно).

Мы оценим модель путем минимизации суммы квадратов ошибок в обоих режимах. Другие подходы, которые учитывают различные отклонения в двух режимах, могут быть лучше, но полезно сохранить их простыми, по крайней мере, на начальном этапе. Следующая функция использует нелинейный оптимизатор для оценки \альфа, \бетаи р. Чтобы обеспечить хорошие начальные значения, мы начинаем оптимизацию с \ Альфа = \ бета = 0и устанавливаем начальное значение ркак среднее из наблюдаемых данных.

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)
      return(1e20)
    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)
    return(sum(e^2))
  }
  fit <- optim(c(0,0,mean(x)),ss,x=x,control=list(maxit=1000))
  if(fit$convergence > 0)
    return(rep(NA,3))
  else
    return(c(alpha=fit$par[1],beta=fit$par[2],r=fit$par[3]))
}

Здесь есть несколько вещей, на которые стоит обратить внимание.

  • Я использовал функцию внутри функции. ssФункция вычисляет сумму квадратов ошибок. Было бы возможно определить это снаружи fitnlts, но, поскольку я не нуждаюсь в этом для каких-либо других целей, лучше определить его локально.
  • Иногда оптимизатор не сходится, а затем fitnltsвозвращает пропущенные значения.
  • Я избегал использования цикла, вычисляя ошибки в обоих режимах. Это намного быстрее, чем альтернативы.
  • Если рза пределами диапазона данных, либо \альфалибо \бетабудет неопознанным. Чтобы избежать этой проблемы, я заставляю рбыть таким, чтобы каждый режим содержал не менее 10% данных. Это произвольное значение, но оно делает оценку более стабильной.
  • Оптимизируемая функция не является непрерывной в рнаправлении. Это приведет к сбою оптимизаторов на основе градиента, но используемый здесь оптимизатор Nelder-Mead относительно устойчив к таким проблемам.
  • Эта конкретная модель может быть более эффективно оценена путем использования условной линейности. Например, мы можем зациклить сетку рзначений (например, в середине между последовательными упорядоченными наблюдениями) и использовать линейную регрессию для других параметров. Но приведенный выше подход носит более общий характер и может быть адаптирован к другим нелинейным моделям временных рядов.

Функции можно использовать следующим образом (с некоторыми примерами параметров):

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

Аналогичный подход может быть использован для других моделей временных рядов.

36
Share: