Identificación de modelos de series temporales con R.



En la fase de identificación aplicaremos distintos procedimientos para presentar modelos tentativos que puedan resumir el comportamiento de la serie temporal.

Incorporamos toda la información disponible sobre la serie:
  • estudios relacionados (la misma serie en otros períodos u otras series que se sospeche tengan comportamiento similar)
  • fenómeno subyacente
  • intervenciones externas
  • variables que modifiquen el comportamiento de nuestra serie
  • intervalo de observación (puede indicar la presencia de periodicidades)
Los pasos para el análisis preliminar de la serie son:

1. Graficar la serie: valores de la serie vs. tiempo. Puede poner de manifiesto:
  • autocorrelación (positiva -persistencia- o negativa -grandes fluctuaciones-)
  • estacionalidad (patrón que se repite aproximadamente transcurridos unos instantes de tiempo)
  • estacionariedad (puede existir no-estacionariedad en la media o en la varianza)
  • tendencias (caso particular de no-estacionariedad en la media)
  • necesidad de transformación de la serie (e.g. Box-Cox)
  • valores extremos (anómalos u outliers)
  • ciclos de período largo
  • intervenciones (cambios en el nivel de la serie por ingerencias externas a la misma)
Realizaremos un ejemplo en R dado por Shumway & Stoffer en su página web (descarcar desde allí el archivo de datos).
#Let's play with the Johnson & Johnson data set. Download jj.dat  to a directory called mydata (or wherever you choose ... the examples below and in the text assume the data are in that directory).
jj = scan("/mydata/jj.dat")        # read the data
jj <- scan("/mydata/jj.dat")       # read the data another way
scan("/mydata/jj.dat") -> jj       # and another
#Now, let's make jj a time series object.
jj = ts(jj, start=1960, frequency=4) #Note that the data are quarterly earnings, hence the frequency=4 statement.

#Now try a plot of the data:
plot(jj, ylab="Earnings per Share", main="J & J")
#Try these and see what happens:
plot(jj, type="o", col="blue", lty="dashed")
plot(diff(log(jj)), main="logged and diffed")

Componentes de la variación: consideraremos los siguientes modelos de descomposición de series temporales.
  • aditivo: Yt=Tt+Ct+St+It. La varianza de ls serie original y la tendencia aumentan con el tiempo. Es decir, St e It son magnitudes que se agregan a Tt y Ct independientemente del valor que tenda.
  • multiplicativo: Yt=Tt*Ct*St+It. El efecto estacional tiene de aumentar al aumentar la tendencia. Es decir, St e It son proporciones de la tendencia Tt.
donde Tt es la tendencia, Ct es la ciclicidad (fluctuaciones en intervalos mayores al año, irregular o al azar), St es la estacionalidad (fluctuaciones periódicas de menso de un año) e It la irregularidad o aleatoriedad o error.

  1. Identificación de la estructura estacionaria: identificar el modelo a ajustar. Lo conseguimos con los pasos anteriores.
  2. Identificación de la estructura no-estacionaria: detectar qué transformaciones hay que aplicar a la serie para conseguir una meda y varianza constantes.
  • Sea xt=mu+yt, podemos eliminar la no-estacionaridad de la tendencia: eliminamos la tendencia al restarle a la serie su media estimada (detrending, yt=xt-mut) o diferenciando la serie (differencing, xt-xt-1). Una ventaja de diferenciar en lugar de eliminar la tendencia es que no estimamos parámetros en la diferenciación, pero como desventaja, la diferenciación no produce una estimación del proceso estocástico yt. Para aplicar diferencias definimos el operador "backshift": Bxt=xt-1, o para cualquierorden d B^dxt=xt-k o diff(d)=(1-B^d).
    Estimar la tendencia: por media móvil, suavizado/filtrado o descomposición. suavizado o alisado (smoothing): por media móvil, kernel de suavizado, regresión ponderado localmente y vecino más cercano, splines.
Ver aquí varios ejemplos de suavizado
  • no-estacionaridad de la varianza: I) transformación logarítimica de la serie: yt=ln(xt). Tiende a eliminar fluctuaciones largas que ocurren sobre porciones de la serie donde los valores subyacentes sin grandes. II) transformación Box-Cox o potencial de la serie: yt= 1) (xt^lambda-1)/lambda si lambda es distinto que cero y 2) ln(xt) si lambda es igual a cero. La transformación Box-Cox se utiliza para mejorar la aproximación a la normalidad o mejorar la linealidad en los valores predichos.
También el análisis de los componentes de variación, el análisis exploratorio de los datos (EDA), el el análisis de regresión y periodograma nos permiten describir la serie.
FILTRADO/SUAVIZADO
#How about filtering/smoothing the Johnson & Johnson series using a two-sided moving average? Let's try this:
fjj(t) = 1/8 jj(t-2) + ¼ jj(t-1) + ¼ jj(t) + ¼ jj(t+1) + 1/8 jj(t+2) #and we'll add a lowess fit for fun.
k = c(.5,1,1,1,.5)            # k is the vector of weights
fjj = filter(jj, sides=2, k)  # ?filter for help [but you knew that already]
plot(jj)
lines(fjj, col="red")         # adds a line to the existing plot
lines(lowess(jj), col="blue", lty="dashed")

DIFERENCIACIÓN
#Let's difference the logged data and call it dljj. Then we'll play with dljj:
dljj = diff(log(jj))        # difference the logged data
plot(dljj)                  # plot it if you haven't already
shapiro.test(dljj)          # test for normality
#Now a histogram and a Q-Q plot, one on top of the other:
par(mfrow=c(2,1))        # set up the graphics
hist(dljj, prob=TRUE, 12)   # histogram
lines(density(dljj))     # smooth it - ?density for details
qqnorm(dljj)             # normal Q-Q plot
qqline(dljj)             # add a line

DESCOMPOSICIÓN
#let's try a structural decomposition of log(jj) = trend + season + error using lowess. Note, this example works only if jj is dimensionless (i.e., you didn't read it in using read.table ... thanks to Jon Moore of the University of Reading, U.K., for pointing this out.)
plot(dog <- stl(log(jj), "per"))
#You can do a little (very little) better using a local seasonal window, plot(dog <- stl(log(jj), s.win=4)), as opposed to the global one used by specifying "per". Type ?stl for details.

REGRESIÓN
#We're going to fit the regression   log(jj)= beta*time + alfa1*Q1 + alfa2*Q2 + alfa3*Q3 + alfa4*Q4 + epsilon where Qi is an indicator of the quarter i = 1,2,3,4. Then we'll inspect the residuals.
Q = factor(rep(1:4,21))     # make (Q)uarter factors [that's repeat 1,2,3,4, 21 times]
trend = time(jj)-1970   # not necessary to "center" time, but the results look nicer
reg = lm(log(jj)~0+trend+Q, na.action=NULL)  # run the regression without an intercept
#-- the na.action statement is to retain time series attributes
summary(reg)
model.matrix(reg) #You can view the model matrix (with the dummy variables) this way:
#Now check out what happened. Look at a plot of the observations and their fitted values:
plot(log(jj), type="o")   # the data in black with little dots
lines(fitted(reg), col=2) # the fitted values in bloody red - or use lines(reg$fitted, col=2)
#and a plot of the residuals and the ACF of the residuals:
par(mfrow=c(2,1))
plot(resid(reg))      # residuals - reg$resid is same as resid(reg)
acf(resid(reg),20)    # acf of the resids. Do those residuals look white? [Ignore the 0-lag correlation, it's always 1.]
#You have to be careful when you regress one time series on lagged components of another using lm(). There is a package called dynlm that makes it easy to fit lagged regressions, and I'll discuss that right after this example. If you use lm(), then what you have to do is "tie" the series together using ts.intersect. If you don't tie the series together, they won't be aligned properly.

2. Función de correlación (simple -FAC- y parcial -FACP-). Sugieren:
  • no-estacionariedad: si el FAC presenta persistencia o decaimiento lento (no exponencial)
  • los cortes/decaimientos en FAC y FACP sugiere la estructura del modelo que debe ajustarse a los datos.
  • matrices de dispersión (scatterplot): observamos las correlaciones a distintos lags.
Conceptos:
  • función media
  • función autocovarianza
  • función de autocorrelación (ACF o FAC)
  • función covarianza cruzada
  • función correlación cruzada (CCF o FCC)
Modelos: ruido blanco, MA(q), AR(p), ARMA()

#Let's check out the correlation structure of dljj using various techniques. First, we'll look at a grid of scatterplots of dljj(t-lag) vs dljj(t) for lag=1,2,...,9.
lag.plot(dljj, 9, do.lines=FALSE)  # why the do.lines=FALSE? ... try leaving it out
#Notice the large positive correlation at lags 4 and 8 and the negative correlations at a few other lags:

#Now let's take a look at the ACF and PACF of dljj:
par(mfrow=c(2,1)) # The power of accurate observation is commonly called cynicism
           #     by those who have not got it. - George Bernard Shaw
acf(dljj, 20)     # ACF to lag 20 - no graph shown... keep reading. Ignore the 0-lag correlation, it's always 1
pacf(dljj, 20)    # PACF to lag 20 - no graph shown... keep reading
#Note that the LAG axis is in terms of frequency, so 1,2,3,4,5 correspond to lags 4,8,12,16,20 because frequency=4 here. If you don't like this type of labeling, you can replace dljj in any of the above by ts(dljj, freq=1); e.g., acf(ts(dljj, freq=1), 20)


3. Modelos preliminares: a partir de los gráficos anteriores, analizamos los componentes de la variación y planteamos todos los modelos posibles que se puedan ajustar a los datos.


Referencias

Comentarios

Publicar un comentario