Fonction de transfert d'intervention ARIMA - Comment visualiser l'effet

11

J'ai une série chronologique mensuelle avec une intervention et je voudrais quantifier l'effet de cette intervention sur le résultat. Je me rends compte que la série est plutôt courte et l'effet n'est pas encore conclu.

Les données

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")

entrez la description de l'image ici

La méthodologie

1) La série pré-intervention (jusqu'en octobre 2013) a été utilisée avec la auto.arimafonction. Le modèle proposé était ARIMA (1,0,0) avec une moyenne non nulle. L'intrigue ACF avait l'air bien.

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) Étant donné le tracé de la série complète, la réponse impulsionnelle a été choisie ci-dessous, avec T = oct 2013,

entrez la description de l'image ici

qui selon cryer et chan peuvent être ajustés comme suit avec la fonction 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

Les résidus de cela sont apparus OK:

entrez la description de l'image ici

L'intrigue des aménagés et réels:

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

entrez la description de l'image ici

Questions

1) Cette méthodologie est-elle correcte pour l'analyse d'intervention?

2) Puis-je regarder l'estimation / SE pour les composantes de la fonction de transfert et dire que l'effet de l'intervention était significatif?

3) Comment visualiser l'effet de la fonction de transfert (le représenter?)

4) Existe-t-il un moyen d'estimer dans quelle mesure l'intervention a augmenté la production après «x» mois? Je suppose que pour cela (et peut-être # 3), je demande comment travailler avec une équation du modèle - s'il s'agissait d'une simple régression linéaire avec des variables muettes (par exemple), je pourrais exécuter des scénarios avec et sans l'intervention et mesurer l'impact - mais je ne sais pas comment travailler ce type de modèle.

AJOUTER

Par requête, voici les résidus des deux paramétrisations.

D'abord de l'ajustement:

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")

entrez la description de l'image ici

Ensuite, à partir de cet ajustement

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")

entrez la description de l'image ici

B_Miner
la source
Serait-il correct si je vous fournissais une solution à l'aide du logiciel SAS?
prévisionniste
Bien sûr, je serais curieux de trouver un meilleur modèle.
B_Miner
D'accord, le modèle est un peu meilleur que celui proposé à l'origine, mais similaire à @javlacalle.
prévisionniste

Réponses:

12

Un modèle AR (1) avec l'intervention définie dans l'équation donnée dans la question peut être ajusté comme indiqué ci-dessous. Remarquez comment l'argument transferest défini; vous avez également besoin d'une variable indicatrice xtransfpour chacune des interventions (le pouls et le changement transitoire):

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

Vous pouvez tester la signification de chaque intervention en examinant la statistique t des coefficients et . Pour plus de commodité, vous pouvez utiliser la fonction .ω 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

Dans ce cas, l'impulsion n'est pas significative au niveau de signification de . Son effet peut déjà être capté par le changement transitoire.5%

L'effet d'intervention peut être quantifié comme suit:

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)

Vous pouvez tracer l'effet de l'intervention comme suit:

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

Effet d'intervention total

L'effet est relativement persistant car est proche de (si était égal à on observerait un décalage de niveau permanent). 1 ω 2 1ω21ω21

Numériquement, ce sont les augmentations estimées quantifiées à chaque instant provoquées par l'intervention en octobre 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

L'intervention augmente la valeur de la variable observée en octobre 2013 d'environ . Dans les périodes suivantes, l'effet persiste mais avec un poids décroissant.75%

Nous pourrions également créer les interventions à la main et les transmettre à stats::arimades régresseurs externes. Les interventions sont une impulsion plus un changement transitoire avec le paramètre et peuvent être construites comme suit.0.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

Les mêmes estimations des coefficients que ci-dessus sont obtenues. Ici, nous avons fixé à . La matrice est le type de variable fictive dont vous pourriez avoir besoin pour essayer différents scénarios. Vous pouvez également définir différentes valeurs pour et comparer son effet. 0,9231 ω 2ω20.9231xregω2

Ces interventions équivalent à une valeur aberrante additive (AO) et à un changement transitoire (TC) définies dans le package tsoutliers. Vous pouvez utiliser ce package pour détecter ces effets comme indiqué dans la réponse de @forecaster ou pour construire les régresseurs utilisés auparavant. Par exemple, dans ce cas:

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

Modifier 1

J'ai vu que l'équation que vous avez donnée peut être réécrite comme suit:

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

et il peut être spécifié comme vous l'avez fait en utilisant transfer=list(c(1, 1)).

Comme indiqué ci-dessous, cette paramétrisation conduit, dans ce cas, à des estimations de paramètres qui impliquent un effet différent par rapport à la paramétrisation précédente. Cela me rappelle l'effet d'une valeur aberrante innovante plutôt qu'une impulsion plus un changement transitoire.

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

Je ne connais pas très bien la notation du package TSAmais je pense que l'effet de l'intervention peut désormais être quantifié comme suit:

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)")

paramétrage de l'effet d'intervention 2

L'effet peut être décrit maintenant comme une forte augmentation en octobre 2013 suivie d'une diminution dans la direction opposée; puis l'effet de l'intervention disparaît rapidement en alternant les effets positifs et négatifs du poids en décomposition.

Cet effet est quelque peu particulier mais peut être possible dans les données réelles. À ce stade, je voudrais examiner le contexte de vos données et les événements qui peuvent avoir affecté les données. Par exemple, y a-t-il eu un changement de politique, une campagne marketing, une découverte, ... qui peut expliquer l'intervention en octobre 2013. Si oui, est-il plus judicieux que cet événement ait un effet sur les données comme décrit précédemment ou comme nous l'avons constaté avec le paramétrage initial?

Selon l'AIC, le modèle initial serait préféré car il est plus bas ( contre ). L'intrigue de la série originale ne suggère pas une correspondance claire avec les changements brusques impliqués dans la mesure de la deuxième variable d'intervention.- 15,4218.9415.42

Sans connaître le contexte des données, je dirais qu'un modèle AR (1) avec un changement transitoire avec le paramètre serait approprié pour modéliser les données et mesurer l'intervention.0.9

Modifier 2

La valeur de détermine la vitesse à laquelle l'effet de l'intervention diminue jusqu'à zéro, c'est donc le paramètre clé du modèle. Nous pouvons l'inspecter en ajustant le modèle pour une plage de valeurs de . Ci-dessous, l'AIC est stocké pour chacun de ces modèles.ω 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 pour différentes valeurs d'oméga

L'AIC le plus bas est trouvé pour (en accord avec la valeur estimée précédemment). Ce paramètre implique un effet relativement persistant mais transitoire. Nous pouvons conclure que l'effet est temporaire car avec des valeurs supérieures à l'AIC augmente (rappelez-vous que dans la limite, , l'intervention devient un décalage de niveau permanent).0,9 ω 2 = 1ω2=0.880.9ω2=1

L'intervention doit être incluse dans les prévisions. L'obtention de prévisions pour des périodes déjà observées est un exercice utile pour évaluer la performance des prévisions. Le code ci-dessous suppose que la série se termine en octobre 2013. Les prévisions sont alors obtenues en incluant l'intervention avec le paramètre .ω2=0.9

On abord le modèle AR (1) avec l'intervention comme régresseur (avec paramètre ):ω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)))

Les prévisions peuvent être obtenues et affichées comme suit:

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")

valeurs observées et prévues

Les premières prévisions correspondent relativement bien aux valeurs observées (ligne pointillée grise). Les prévisions restantes montrent comment la série poursuivra son chemin vers la moyenne d'origine. Les intervalles de confiance sont néanmoins importants, reflétant l'incertitude. Nous devons donc être prudents et réviser le modèle à mesure que de nouvelles données sont enregistrées.

95% intervalles de confiance de peuvent être ajoutés au tracé précédent comme suit:

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
la source
C'est super, merci! J'ai eu quelques suivis si cela ne vous dérange pas. 1) Le processus que j'ai suivi est-il correct? 2) Estimeriez-vous que l'ajustement du modèle est «suffisamment bon» pour utiliser les estimations afin de quantifier l'effet de l'intervention? 3) Ne devrais-je pas pouvoir utiliser ma paramétrisation, c'est-à-dire transfer = list (c (1,1)) comme équivalent et obtenir des résultats assez proches? L'exemple que je
suivais
@B_Miner Vous avez raison, j'ai modifié ma réponse.
javlacalle
Je suis d'accord avec vous que des deux modèles, la première paramétrisation (peut-être avec l'impulsion non significative supprimée) serait la meilleure solution. Pourquoi les deux paramétrisations ne donnent pas le même modèle, alors que je le pense, est un mystère. Je vais envoyer un e-mail au développeur du package (qui a également écrit le livre qui mentionne leur équivalence).
B_Miner
Les données sont le nombre de certificats de dépôts ouverts par mois. L'intervention a été une augmentation du taux d'intérêt moyen, qui a augmenté à partir du 13 octobre. Le niveau du taux d'intérêt est resté relativement constant depuis le 13 octobre. Il m'a semblé qu'après la flambée, la demande pour le produit a commencé à diminuer - Je ne sais pas s'il reviendra à la moyenne précédente ou s'installera à un niveau élevé (par rapport au précédent).
B_Miner
B_miner, sur la base des données que nous ne pouvons pas vraiment conclure, si la demande se stabilise à une nouvelle moyenne.
prévisionniste
4

Parfois, moins c'est plus. Avec 30 observations en main, j'ai soumis les données à AUTOBOX, un logiciel que j'ai aidé à développer. Je soumets l'analyse suivante dans l'espoir d'obtenir la récompense +200 (je plaisante!). J'ai tracé les valeurs réelles et nettoyées en suggérant visuellement l'impact de "l'activité récente". entrez la description de l'image ici. Le modèle développé automatiquement est présenté ici. entrez la description de l'image iciet ici entrez la description de l'image ici. Les résidus de cette série assez simple à décalage de niveau sont présentés ici entrez la description de l'image ici. Les statistiques du modèle sont ici entrez la description de l'image ici. En résumé, il y avait une intervention qui pourrait être identifiée empiriquement rendant un processus ARIMA; deux impulsions et 1 décalage de niveau entrez la description de l'image ici. Le graphique Réel / Ajustement et Prévisions met davantage en évidence l'analyse.entrez la description de l'image ici

Pour ma part, je voudrais voir le tracé des résidus des modèles précédemment spécifiés et à mon avis potentiellement trop spécifiés.

IrishStat
la source
Je ne connais pas Autobox, mais la partie bruit du modèle est-elle la même que celle que j'avais à l'origine: une moyenne non nulle et un AR (1)?
B_Miner
Est-ce que cette sortie indique que la seule "intervention" dans la période du 13 octobre à la période actuelle est une seule impulsion pour le 13 octobre, puis la série revient à son niveau moyen normal?
B_Miner
J'ai ajouté les résidus des deux paramétrisations. À mes yeux, il semble que le premier que j'ai énuméré (celui initialement adapté par javlacalle) soit meilleur. Se mettre d'accord?
B_Miner
1) La partie bruit est un AR (1) avec une moyenne non nulle
IrishStat
1) La partie bruit est un AR (1) avec une moyenne non nulle; 2) Il y a 2 interventions période 22 et période 3 et après le 13 octobre il revient à un nouveau niveau qui a commencé au 13 septembre; 3) Étant donné le choix entre les deux que vous avez mentionnés, je suis d'accord MAIS je préfère le modèle AUTOBOX pour sa simplicité et son efficacité. Vous pouvez en savoir plus sur AUTOBOX sur autobox.com/cms
IrishStat
3

R

Voici le code:

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

Ci-dessous l'estimation, il y a eu une augmentation de ~ 2356,3 unités en octobre 2013 avec une erreur standard de ~ 481,8 et a un effet de décomposition par la suite. La fonction a identifié automatiquement AR (1). J'ai dû faire quelques itérations et faire la différence saisonnière et non saisonnière à 0, ce qui se reflète dans la méthode args.ts dans la fonction 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

Ci-dessous est l'intrigue, tsoutlier est le seul package que je connaisse qui puisse imprimer ces modifications temporaires dans un graphique.

entrez la description de l'image ici

Nous espérons que cette analyse a répondu à vos 2, 3 et 4 questions, mais en utilisant une méthdéologie différente. Surtout le graphique et les coefficients ont fourni l'effet de cette intervention et ce qui se serait passé si vous n'aviez pas eu cette intervention.

En espérant également que quelqu'un d'autre puisse reproduire cette intrigue / analyse en utilisant la modélisation de la fonction de transfert dans R. Je ne sais pas si cela pourrait être fait dans R, peut-être que quelqu'un d'autre peut me vérifier à ce sujet.

prévisionniste
la source
Merci. Oui, cette intrigue est ce que je voudrais faire à partir du modèle arimax - regardez avec et sans l'intervention (et soustrayez). Je pense que la fonction de filtre dans R peut être utilisée pour générer la valeur de la fonction de transfert pour chaque mois (et ensuite simplement la tracer pour la visualiser), mais je n'arrive pas à comprendre comment le faire pour une fonction d'intervention par impulsion arbitraire.
B_Miner