Измерение характеристик временных рядов

Несколько лет назад я работал над проектом, в котором мы измеряли различные характеристики временного ряда и использовали информацию, чтобы определить, какой метод прогнозирования применить или как сгруппировать временные ряды в значимые группы. Два основных документа, которые вышли из этого проекта:

С тех пор у меня было много запросов на код, который один из моих соавторов отправил по электронной почте всем, кто спросил. Но чтобы сделать это проще, мы подумали, что будет полезно, если я выложу здесь обновленный код. Это не то же самое, что R-код, который мы использовали в статье, поскольку я улучшил его несколькими способами (поэтому он даст разные результаты). Если вы просто хотите код, перейдите к нижней части поста.

Нахождение периода данных

Обычно при работе с временными рядами мы знаем период данных (например, если наблюдения ежемесячные, период равен 12). Но в этом проекте некоторые из наших данных были неизвестного периода, и мы хотели, чтобы метод автоматически определил соответствующий период. Метод, который мы использовали, был основан на локальных пиках и впадинах в АКФ . Но с тех пор я разработал лучший подход ( предложенный на перекрестной валидации .com ) с использованием оценки спектральной плотности:

find.freq <- function(x)
{
  n <- length(x)
  spec <- spec.ar(c(na.contiguous(x)),plot=FALSE)
  if(max(spec$spec)>10) # Arbitrary threshold chosen by trial and error.
  {
    period <- round(1/spec$freq[which.max(spec$spec)])
    if(period==Inf) # Find next local maximum
    {
      j <- which(diff(spec$spec)>0)
      if(length(j)>0)
      {
        nextmax <- j[1] + which.max(spec$spec[j[1]:500])
        if(nextmax <= length(spec$freq))
          period <- round(1/spec$freq[nextmax])
        else
          period <- 1
      }
      else
        period <- 1
    }
  }
  else
    period <- 1
 
  return(period)
}

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

Разложение данных на трендовую и сезонную составляющие

Нам нужно было измерить силу тренда и силу сезонности, и для этого мы разбили данные на тренды, сезонность и ошибки.

Поскольку не все данные можно было аддитивно разложить, нам сначала нужно было применить автоматическое преобразование Бокса-Кокса. Мы опробовали диапазон параметров Бокса-Кокса на сетке и выбрали тот, который дает наиболее нормальные ошибки. Это работало хорошо, но с тех пор я нашел несколько статей, которые предоставляют довольно хорошие автоматизированные алгоритмы Бокса-Кокса, которые я реализовал в пакете прогноза. Так что в этом коде вместо этого используется метод Герреро (1993).

Для сезонных временных рядов мы разложили преобразованные данные, используя разложение stl с периодической сезонностью.

Для не- сезонных временных рядов, мы оценили тенденцию преобразованных данных с использованием наказываться регрессии сплайнов с помощью пакета mgcv .

decomp <- function(x,transform=TRUE)
{
  require(forecast)
  # Transform series
  if(transform & min(x,na.rm=TRUE) >= 0)
  {
    lambda <- BoxCox.lambda(na.contiguous(x))
    x <- BoxCox(x,lambda)
  }
  else
  {
    lambda <- NULL
    transform <- FALSE
  }
  # Seasonal data
  if(frequency(x)>1)
  {
    x.stl <- stl(x,s.window="periodic",na.action=na.contiguous)
    trend <- x.stl$time.series[,2]
    season <- x.stl$time.series[,1]
    remainder <- x - trend - season
  }
  else #Nonseasonal data
  {
    require(mgcv)
    tt <- 1:length(x)
    trend <- rep(NA,length(x))
    trend[!is.na(x)] <- fitted(gam(x ~ s(tt)))
    season <- NULL
    remainder <- x - trend
  }
  return(list(x=x,trend=trend,season=season,remainder=remainder,
    transform=transform,lambda=lambda))
}

Расставляем все по шкале [0,1]

Мы хотели измерить ряд характеристик, таких как сила сезонности, сила тренда, уровень нелинейности, асимметрия, эксцесс, серийная корреляция, самоподобие, уровень хаотичности (это слово?) И периодичность данные. Но мы хотели все это в одном масштабе, что означало отображение естественного диапазона каждой меры на [0,1]. Следующие две функции были использованы для этого.

# f1 maps [0,infinity) to [0,1]
f1 <- function(x,a,b)
{
  eax <- exp(a*x)
  if (eax == Inf)
    f1eax <- 1
  else
    f1eax <- (eax-1)/(eax+b)
  return(f1eax)
}
 
# f2 maps [0,1] onto [0,1]
f2 <- function(x,a,b)
{
  eax <- exp(a*x)
  ea <- exp(a)
  return((eax-1)/(eax+b)*(ea+b)/(ea-1))
}

Значения $a$ и $б$ в каждой функции были выбраны таким образом, чтобы показатель имел 90-й процентиль 0,10, когда данные были стандартно нормальными, и значение 0,9 с использованием хорошо известных эталонных временных рядов.

Расчет мер

Теперь мы готовы рассчитать показатели по исходным данным, а также по скорректированным данным (после удаления тренда и сезонности).

measures <- function(x)
{
  require(forecast)
 
  N <- length(x)
  freq <- find.freq(x)
  fx <- c(frequency=(exp((freq-1)/50)-1)/(1+exp((freq-1)/50)))
  x <- ts(x,f=freq)
 
  # Decomposition
  decomp.x <- decomp(x)
 
  # Adjust data
  if(freq > 1)
    fits <- decomp.x$trend + decomp.x$season
  else # Nonseasonal data
    fits <- decomp.x$trend
  adj.x <- decomp.x$x - fits + mean(decomp.x$trend, na.rm=TRUE)
 
  # Backtransformation of adjusted data
  if(decomp.x$transform)
    tadj.x <- InvBoxCox(adj.x,decomp.x$lambda)
  else
    tadj.x <- adj.x
 
  # Trend and seasonal measures
  v.adj <- var(adj.x, na.rm=TRUE)
  if(freq > 1)
  {
    detrend <- decomp.x$x - decomp.x$trend
    deseason <- decomp.x$x - decomp.x$season
    trend <- ifelse(var(deseason,na.rm=TRUE) < 1e-10, 0, 
      max(0,min(1,1-v.adj/var(deseason,na.rm=TRUE))))
    season <- ifelse(var(detrend,na.rm=TRUE) < 1e-10, 0,
      max(0,min(1,1-v.adj/var(detrend,na.rm=TRUE))))
  }
  else #Nonseasonal data
  {
    trend <- ifelse(var(decomp.x$x,na.rm=TRUE) < 1e-10, 0,
      max(0,min(1,1-v.adj/var(decomp.x$x,na.rm=TRUE))))
    season <- 0
  }
 
  m <- c(fx,trend,season)
 
  # Measures on original data
  xbar <- mean(x,na.rm=TRUE)
  s <- sd(x,na.rm=TRUE)
 
  # Serial correlation
  Q <- Box.test(x,lag=10)$statistic/(N*10)
  fQ <- f2(Q,7.53,0.103)
 
  # Nonlinearity
  p <- terasvirta.test(na.contiguous(x))$statistic
  fp <- f1(p,0.069,2.304)
 
  # Skewness
  sk <- abs(mean((x-xbar)^3,na.rm=TRUE)/s^3)
  fs <- f1(sk,1.510,5.993)
 
  # Kurtosis
  k <- mean((x-xbar)^4,na.rm=TRUE)/s^4
  fk <- f1(k,2.273,11567)
 
  # Hurst=d+0.5 where d is fractional difference.
  H <- fracdiff(na.contiguous(x),0,0)$d + 0.5
 
  # Lyapunov Exponent
  if(freq > N-10)
      stop("Insufficient data")
  Ly <- numeric(N-freq)
  for(i in 1:(N-freq))
  {
    idx <- order(abs(x[i] - x))
    idx <- idx[idx < (N-freq)]
    j <- idx[2]
    Ly[i] <- log(abs((x[i+freq] - x[j+freq])/(x[i]-x[j])))/freq
    if(is.na(Ly[i]) | Ly[i]==Inf | Ly[i]==-Inf)
      Ly[i] <- NA
  }
  Lyap <- mean(Ly,na.rm=TRUE)
  fLyap <- exp(Lyap)/(1+exp(Lyap))
 
  m <- c(m,fQ,fp,fs,fk,H,fLyap)
 
  # Measures on adjusted data
  xbar <- mean(tadj.x, na.rm=TRUE)
  s <- sd(tadj.x, na.rm=TRUE)
 
  # Serial
  Q <- Box.test(adj.x,lag=10)$statistic/(N*10)
  fQ <- f2(Q,7.53,0.103)
 
  # Nonlinearity
  p <- terasvirta.test(na.contiguous(adj.x))$statistic
  fp <- f1(p,0.069,2.304)
 
  # Skewness
  sk <- abs(mean((tadj.x-xbar)^3,na.rm=TRUE)/s^3)
  fs <- f1(sk,1.510,5.993)
 
  # Kurtosis
  k <- mean((tadj.x-xbar)^4,na.rm=TRUE)/s^4
  fk <- f1(k,2.273,11567)
 
  m <- c(m,fQ,fp,fs,fk)
  names(m) <- c("frequency", "trend","seasonal",
    "autocorrelation","non-linear","skewness","kurtosis",
    "Hurst","Lyapunov",
    "dc autocorrelation","dc non-linear","dc skewness","dc kurtosis")
 
  return(m)
}

Вот краткий пример применения ежемесячной добычи газа в Австралии:

library(forecast)
measures(gas)
     frequency              trend           seasonal    autocorrelation 
        0.1096             0.9989             0.9337             0.9985 
    non-linear           skewness           kurtosis              Hurst 
        0.4947             0.1282             0.0055             0.9996 
      Lyapunov dc autocorrelation      dc non-linear        dc skewness 
        0.5662             0.1140             0.0538             0.1743 
   dc kurtosis 
        0.9992

Функция далека от совершенства, и нетрудно найти примеры, где она не работает. Например, он не работает с множественной сезонностью — попробуйте измерить (Тейлор) и проверить сезонность. Кроме того, я не уверен, что эксцесс дает здесь что-то полезное или что измерение асимметрии выполняется наилучшим образом. Но это действительно было доказательством концепции, поэтому мы предоставим другим возможность пересматривать и улучшать код.

В наших работах, мы приняли меры , полученные с использованием R, и производятся самоорганизующиеся карты с использованием Viscovery . Теперь для этого есть пакет som в R, так что может быть возможно интегрировать и этот шаг в R. Dendogram был создан в MATLAB, хотя это может теперь быть сделано в R , используя ggdendro пакет , например.

Загрузите код в одном файле.

Измерение характеристик временных рядов

Нахождение периода данных

Разложение данных на трендовую и сезонную составляющие

Расставляем все по шкале [0,1]

Расчет мер

Категории

Последние статьи

Рефакторинг Hudson God Class

Альтернативы синтаксиса Java лямбда

Morphia и MongoDB: развивающиеся структуры документов

OpenShift Express: развертывание приложения Java EE (с поддержкой AS7)

Интеграция jqGrid, REST, AJAX и Spring MVC