Функция передачи вмешательства ARIMA - как визуализировать эффект

11

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

Данные

cds <- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 3362L,
                   2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L,
                   2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L,
                   4523L, 4186L, 4070L, 4000L, 3498L),
                 .Dim=c(29L, 1L),
                 .Dimnames=list(NULL, "CD"),
                 .Tsp=c(2012, 2014.33333333333, 12), class="ts")

введите описание изображения здесь

Методология

1) С этой auto.arimaфункцией была использована серия до вмешательства (до октября 2013 года) . Предложенная модель была ARIMA (1,0,0) с ненулевым средним. Сюжет ACF выглядел хорошо.

pre <- window(cds, start=c(2012, 01), end=c(2013, 09))

mod.pre <- auto.arima(log(pre))

# Coefficients:
#          ar1  intercept
#       0.5821     7.9652
# s.e.  0.1763     0.0810
# 
# sigma^2 estimated as 0.02709:  log likelihood=7.89
# AIC=-9.77   AICc=-8.36   BIC=-6.64

2) Учитывая график полной серии, импульсный отклик был выбран ниже, с T = октябрь 2013,

введите описание изображения здесь

который в соответствии с cryer и chan может быть подобран следующим образом с помощью функции arimax:

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1)))
mod.arimax

# Series: log(cds) 
# ARIMA(1,0,0) with non-zero mean 
# 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# 
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

Остатки от этого появились ОК:

введите описание изображения здесь

Сюжет обустроен и актуален:

plot(fitted(mod.arimax), col="red", type="b")
lines(window(log(cds), start=c(2012, 02)), type="b")

введите описание изображения здесь

Вопросы

1) Является ли эта методология правильной для анализа вмешательства?

2) Могу ли я посмотреть на оценку / SE для компонентов передаточной функции и сказать, что эффект вмешательства был значительным?

3) Как можно визуализировать эффект передаточной функции (построить ее?)

4) Есть ли способ оценить, насколько вмешательство увеличило результат после «х» месяцев? Я предполагаю для этого (и, возможно, № 3), я спрашиваю, как работать с уравнением модели - если бы это была простая линейная регрессия с фиктивными переменными (например), я мог бы запустить сценарии с вмешательством и без вмешательства и измерить воздействие - но я просто не уверен, как работать с этим типом модели.

ДОБАВИТЬ

По запросу, здесь остатки от двух параметризаций.

Первый из подгонки:

fit <- arimax(log(cds), order=c(1, 0, 0),
              xtransf=
              data.frame(Oct13a=1 * (seq_along(cds) == 22),
                         Oct13b=1 * (seq_along(cds) == 22)),
              transfer=list(c(0, 0), c(1, 0)))

plot(resid(fit), type="b")

введите описание изображения здесь

Тогда из этого подходит

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1))) 

mod.arimax
plot(resid(mod.arimax), type="b")

введите описание изображения здесь

B_Miner
источник
Было бы хорошо, если бы я предоставил вам решение с использованием программного обеспечения SAS?
синоптик
Конечно, мне было бы интересно, если бы вы придумали лучшую модель.
B_Miner
Хорошо, модель немного лучше, чем предполагалось изначально, но похожа на @javlacalle.
синоптик

Ответы:

12

Модель AR (1) с вмешательством, определенным в уравнении, приведенном в вопросе, можно подобрать, как показано ниже. Обратите внимание, как определяется аргумент transfer; Вам также нужна одна индикаторная переменная xtransfдля каждого из вмешательств (импульс и временное изменение):

require(TSA)
cds <- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 3362L,
                   2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L,
                   2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L,
                   4523L, 4186L, 4070L, 4000L, 3498L),
                 .Dim = c(29L, 1L),
                 .Dimnames = list(NULL, "CD"),
                 .Tsp = c(2012, 2014.33333333333, 12),
                 class = "ts")

fit <- arimax(log(cds), order = c(1, 0, 0), 
              xtransf = data.frame(Oct13a = 1 * (seq_along(cds) == 22), 
                                   Oct13b = 1 * (seq_along(cds) == 22)),
              transfer = list(c(0, 0), c(1, 0)))
fit
# Coefficients:
#          ar1  intercept  Oct13a-MA0  Oct13b-AR1  Oct13b-MA0
#       0.5599     7.9643      0.1251      0.9231      0.4332
# s.e.  0.1563     0.0684      0.1911      0.1146      0.2168
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -18.94

Вы можете проверить значимость каждого вмешательства, посмотрев на t-статистику коэффициентов и . Для удобства вы можете использовать функцию .ω 1ω0ω1coeftest

require(lmtest)
coeftest(fit)
#            Estimate Std. Error  z value  Pr(>|z|)    
# ar1        0.559855   0.156334   3.5811 0.0003421 ***
# intercept  7.964324   0.068369 116.4896 < 2.2e-16 ***
# Oct13a-MA0 0.125059   0.191067   0.6545 0.5127720    
# Oct13b-AR1 0.923112   0.114581   8.0564 7.858e-16 ***
# Oct13b-MA0 0.433213   0.216835   1.9979 0.0457281 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

В этом случае пульс незначителен на уровне значимости . Его эффект может быть уже уловлен временным изменением.5%

Эффект вмешательства может быть определен следующим образом:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(
  intv.effect * 0.1251 + 
  filter(intv.effect, filter = 0.9231, method = "rec", sides = 1) * 0.4332)
intv.effect <- exp(intv.effect)
tsp(intv.effect) <- tsp(cds)

Вы можете построить эффект вмешательства следующим образом:

plot(100 * (intv.effect - 1), type = "h", main = "Total intervention effect")

Общий эффект вмешательства

1ω21ω21

Численно, это предполагаемые увеличения, количественно определенные в каждый момент времени, вызванные вмешательством в октябре 2013 года:

window(100 * (intv.effect - 1), start = c(2013, 10))
#           Jan      Feb      Mar      Apr      May Jun Jul Aug Sep      Oct
# 2013                                                              74.76989
# 2014 40.60004 36.96366 33.69046 30.73844 28.07132                         
#           Nov      Dec
# 2013 49.16560 44.64838

75%

stats::arima0.9231

xreg <- cbind(
  I1 = 1 * (seq_along(cds) == 22), 
  I2 = filter(1 * (seq_along(cds) == 22), filter = 0.9231, method = "rec", 
              sides = 1))
arima(log(cds), order = c(1, 0, 0), xreg = xreg)
# Coefficients:
#          ar1  intercept      I1      I2
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -20.94

ω20.9231xregω2

Эти вмешательства эквивалентны аддитивному выбросу (AO) и временному изменению (TC), определенному в пакете tsoutliers. Вы можете использовать этот пакет для обнаружения этих эффектов, как показано в ответе @forecaster, или для построения регрессоров, использовавшихся ранее. Например, в этом случае:

require(tsoutliers)
mo <- outliers(c("AO", "TC"), c(22, 22))
oe <- outliers.effects(mo, length(cds), delta = 0.9231)
arima(log(cds), order = c(1, 0, 0), xreg = oe)
# Coefficients:
#          ar1  intercept    AO22    TC22
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood=14.47
# AIC=-20.94   AICc=-18.33   BIC=-14.1

Редактировать 1

Я видел, что уравнение, которое вы дали, можно переписать так:

(ω0+ω1)ω0ω2B1ω2BPt

и это может быть указано, как вы использовали transfer=list(c(1, 1)).

Как показано ниже, эта параметризация приводит, в этом случае, к оценкам параметров, которые включают другой эффект по сравнению с предыдущей параметризацией. Это напоминает мне эффект инновационного выброса, а не импульса плюс временное изменение.

fit2 <- arimax(log(cds), order=c(1, 0, 0), include.mean = TRUE, 
  xtransf=data.frame(Oct13 = 1 * (seq(cds) == 22)), transfer = list(c(1, 1)))
fit2
# ARIMA(1,0,0) with non-zero mean 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

Я не очень знаком с обозначением пакета, TSAно я думаю, что эффект от вмешательства теперь может быть количественно определен следующим образом:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(intv.effect * 0.4261 + 
  filter(intv.effect, filter = -0.4429, method = "rec", sides = 1) * 0.3567)
tsp(intv.effect) <- tsp(cds)
window(100 * (exp(intv.effect) - 1), start = c(2013, 10))
#              Jan         Feb         Mar         Apr         May Jun Jul Aug
# 2014  -3.0514633   1.3820052  -0.6060551   0.2696013  -0.1191747            
#      Sep         Oct         Nov         Dec
# 2013     118.7588947 -14.6135216   7.2476455

plot(100 * (exp(intv.effect) - 1), type = "h", 
  main = "Intervention effect (parameterization 2)")

параметризация эффекта вмешательства 2

Эффект теперь можно описать как резкое увеличение в октябре 2013 года с последующим снижением в противоположном направлении; затем эффект вмешательства быстро исчезает, чередуя положительные и отрицательные эффекты от снижения веса.

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

18.9415.42

0.9

Редактировать 2

ω2ω2

omegas <- seq(0.5, 1, by = 0.01)
aics <- rep(NA, length(omegas))
for (i in seq(along = omegas)) {
  tc <- filter(1 * (seq_along(cds) == 22), filter = omegas[i], method = "rec", 
               sides = 1)
  tc <- ts(tc, start = start(cds), frequency = frequency(cds))
  fit <- arima(log(cds), order = c(1, 0, 0), xreg = tc)
  aics[i] <- AIC(fit)
}
omegas[which.min(aics)]
# [1] 0.88

plot(omegas, aics, main = "AIC for different values of the TC parameter")

AIC для разных значений омега

ω2=0.880.9ω2=1

ω2=0.9

ω2=0.9

tc <- filter(1 * (seq.int(length(cds) + 12) == 22), filter = 0.9, method = "rec", 
             sides = 1)
tc <- ts(tc, start = start(cds), frequency = frequency(cds))
fit <- arima(window(log(cds), end = c(2013, 10)), order = c(1, 0, 0), 
             xreg = window(tc, end = c(2013, 10)))

Прогнозы могут быть получены и отображены следующим образом:

p <- predict(fit, n.ahead = 19, newxreg = window(tc, start = c(2013, 11)))

plot(cbind(window(cds, end = c(2013, 10)), exp(p$pred)), plot.type = "single", 
     ylab = "", type = "n")
lines(window(cds, end = c(2013, 10)), type = "b")
lines(window(cds, start = c(2013, 10)), col = "gray", lty = 2, type = "b")
lines(exp(p$pred), type = "b", col = "blue")
legend("topleft",
       legend = c("observed before the intervention",
           "observed after the intervention", "forecasts"),
       lty = rep(1, 3), col = c("black", "gray", "blue"), bty = "n")

наблюдаемые и прогнозируемые значения

Первые прогнозы относительно хорошо соответствуют наблюдаемым значениям (серая пунктирная линия). Остальные прогнозы показывают, как ряд будет продолжать путь к исходному среднему значению. Доверительные интервалы, тем не менее, велики, что отражает неопределенность. Поэтому следует проявлять осторожность и пересматривать модель по мере записи новых данных.

95%

lines(exp(p$pred + 1.96 * p$se), lty = 2, col = "red")
lines(exp(p$pred - 1.96 * p$se), lty = 2, col = "red")
javlacalle
источник
Это здорово, спасибо! У меня была пара продолжений, если вы не возражаете. 1) Правильный ли процесс, которому я следовал? 2) Считаете ли вы соответствие модели «достаточно хорошим», чтобы использовать оценки для количественной оценки эффекта вмешательства? 3) Должен ли я не быть в состоянии использовать свою параметризацию, то есть Transfer = List (c (1,1)) в качестве эквивалента и получить довольно близкие результаты? Пример, которому я следовал из учебника, предполагал, что я должен быть в состоянии, но в этом примере результаты не близки ...
B_Miner
@B_Miner Вы правы, я отредактировал свой ответ.
Javlacalle
Я согласен с вами в том, что из двух моделей первая параметризация (возможно, с удаленным незначительным импульсом) будет наиболее подходящей. Почему две параметризации не дают одну и ту же модель, когда я полагаю, что это должно быть, остается загадкой. Я отправлю электронное письмо разработчику пакета (который также написал книгу, в которой упоминается их эквивалентность).
B_Miner
Данные - количество депозитных сертификатов, открытых в месяц. Интервенция заключалась в увеличении средней процентной ставки, которая начала расти 13 октября. Уровень процентной ставки оставался относительно постоянным с 13 октября. Мне показалось, что после пика спрос на продукт начал снижаться - Я не уверен, вернется ли он к предыдущему среднему значению или установится на каком-то повышенном (от предыдущего) уровне.
B_Miner
B_miner, основываясь на данных, которые мы не можем на самом деле сделать вывод, если спрос снизится до нового среднего значения.
синоптик
4

Иногда меньше значит больше. Имея 30 наблюдений, я отправил данные в AUTOBOX, программное обеспечение, которое помогло мне разработать. Я представляю следующий анализ в надежде получить награду +200 (шучу!). Я составил фактические и очищенные значения, наглядно демонстрируя влияние «недавних действий». введите описание изображения здесь, Модель, которая была разработана автоматически, показана здесь. введите описание изображения здесьи здесь введите описание изображения здесь. Остатки от этого довольно простого смещенного уровня представлены здесь введите описание изображения здесь. Модель статистики здесь введите описание изображения здесь. Таким образом, были вмешательства, которые можно было бы определить эмпирически, представляя процесс ARIMA; два импульса и 1 сдвиг уровня введите описание изображения здесь. График Actual / Fit и Forecast дополнительно выделяет анализ.введите описание изображения здесь

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

IrishStat
источник
Я не знаком с Autobox, но шумовая часть модели такая же, как у меня изначально: ненулевое среднее и AR (1)?
B_Miner
Говорит ли этот вывод, что единственное «вмешательство» в период с 13 октября по текущий период времени - это одиночный импульс для 13 октября, а затем серия возвращается к своему нормальному среднему уровню?
B_Miner
Я добавил остатки от обеих параметризаций. На мой взгляд, кажется, что первый, который я перечислил (тот, который изначально был у javlacalle), лучше. Согласен?
B_Miner
1)
Шумовая
1) шумовой частью является AR (1) с ненулевым средним; 2) Есть 2 вмешательства, период 22 и период 3, и после 13 октября он возвращается на новый уровень, который начался 13 сентября; 3) Учитывая выбор между двумя, которые вы упомянули, я согласен, НО я предпочитаю модель AUTOBOX за ее простоту и эффективность. Вы можете узнать больше об AUTOBOX от autobox.com/cms
IrishStat
3

R

Ниже приведен код:

cds<- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L, 
                  3362L, 2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L, 
                  2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L, 4523L, 
                  4186L, 4070L, 4000L, 3498L), .Dim = c(29L, 1L), .Dimnames = list(
                    NULL, "CD"), .Tsp = c(2012, 2014.33333333333, 12), class = "ts")
arimatr <- tsoutliers::tso(cds,args.tsmethod=list(d=0,D=0))
plot(arimatr)
arimatr

Ниже приведена оценка, что в октябре 2013 года наблюдалось увеличение на ~ 2356,3 единиц со стандартной ошибкой ~ 481,8, и в последующем он имеет затухающий эффект. Функция автоматически идентифицирует AR (1). Мне пришлось сделать пару итераций и сделать как сезонное, так и несезонное различие до 0, что отражено в args.tsmethod в функции tso.

Series: cds 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1  intercept       TC22
      0.5969  3034.6560  2356.2914
s.e.  0.1495   206.5202   481.7981

sigma^2 estimated as 209494:  log likelihood=-219.03
AIC=446.06   AICc=447.73   BIC=451.53

Outliers:
  type ind    time coefhat tstat
1   TC  22 2013:10    2356 4.891

Ниже приведен график, tsoutlier - единственный известный мне пакет, который может печатать временные изменения в графике.

введите описание изображения здесь

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

Также надеясь, что кто-то еще сможет воспроизвести этот график / анализ, используя моделирование передаточной функции в R. Я не уверен, что это может быть сделано в R, может быть кто-то еще может проверить меня на этот счет.

предсказатель
источник
Спасибо. Да, этот график - то, что я хотел бы сделать из модели arimax - смотреть с вмешательством и без него (и вычитать). Я думаю, что функцию фильтра в R можно использовать для генерации значения передаточной функции для каждого месяца (а затем просто построить его для визуализации), но я не могу понять, как это сделать для произвольной импульсной интервенционной функции.
B_Miner