Principios y Aplicaciones

Mate científicoMate científico "At least once in a lifetime it is convenient to put everything to discussion" (Descartes)

Readmore

useRLengua franca Usted no puede, no debe dejar de tenerlo...

Readmore

SuperCarlitosSuperCarlitos Y pa mi que la evolución tuvo que empesá má o meno así...

Readmore

Libros para descargar sobre Análisis de Series Temporales

martes, 11 de enero de 2011

0
Análisis de series de temperatura global en R.

### DATOS: debemos descargar los datos desde http://www.stat.pitt.edu/stoffer/tsa2/index.html y guardarlos en el directorio de trabajo de R. 
#gtemp - Global mean land-ocean temperature deviations (from 1951-1980 average), measured in degrees centigrade, for the years 1880-2009;data taken from http://data.giss.nasa.gov/gistemp/graphs/
#gtemp2 - Similar to gtemp but the data are based only on surface air temperature data obtained from meteorological stations.

############################# 2º edición del libro ################
#Example 1.2
globtemp = ts(scan("globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900) # use data from 1990
plot(gtemp, type="o", ylab="Global Temperature Deviations")

#Example 2.1
fit = lm(gtemp~time(gtemp), na.action=NULL) # regress gtemp on time
#-- note: na.action=NULL is used to retain time series attributes - ?lm for info
summary(fit) # view the results
plot(gtemp, type="o", ylab="Global Temperature Deviation") # plot the series
abline(fit) # add regression line to the plot

#Examples 2.3 and 2.4
fit = lm(gtemp~time(gtemp), na.action=NULL) # regress gtemp on time

x11(); par(mfrow=c(3,1))
plot(gtemp, type="o", ylab="Global Temperature Deviation"); abline(fit) # add regression line to the plot
plot(resid(fit), type="o", main="detrended")
plot(diff(gtemp), type="o", main="first difference") # this and below will do example 2.4


x11() ; par(mfrow=c(3,1)) # plot ACFs
acf(gtemp, 48)
acf(resid(fit), 48)
acf(diff(gtemp), 48)


0
Comandos importantes. (Stoffer)

ISSUE 5: When is lead, lag, and why is yesterday perfectly correlated with today?

You have to be careful when working with lagged components of a time series. Note that lag(x,1) is a FORWARD shift and lag(x,-1) is a BACKWARD shift. Try a small example:

x=ts(1:5)
cbind(x,lag(x,1),lag(x,-1))

Time Series:
Start = 0
End = 6
Frequency = 1
x lag(x, 1) lag(x, -1)
0 NA 1 NA
1 1 2 NA
2 2 3 1
3 3 4 2 <-- here, x is 3, lag(x,1) is 4, lag(x,-1) is 2
4 4 5 3
5 5 NA 4
6 NA NA 5

Note that the default of the command lag(x) is lag(x,1). So, in R, if you have a series x(t), then
y(t) = lag{x(t)} = x(t+1), and NOT x(t-1).
This seems awkward, and it's not typical of other programs. But, that's the way Splus does it, so why not R? As long as you know the convention, you'll be ok ...

... well, then I started wondering how this plays out in other things. So, I started playing around with some commands. In what you'll see next, I'm using two simultaneously measured series presented in the text called soi and rec... it doesn't matter what they are for this demonstration. First, I entered the command

acf(cbind(soi, rec))

and I got:

ccf0

Before you scroll down, try to figure out what the graphs are giving you (in particular, the off-diagonal plots ... and yes they're CCFs, but what's the lead-lag relationship in each plot???) ...
.
.
.
.
.
.
.
... here you go:

ccf

The jpg is messy, but you'll get the point... the writing is mine. When you see something like 'rec "leads" here', it means rec comes in time before soi, and so on. Anyway, to be consistent, shouldn't the graph in the 2nd row, 1st column be corr{rec(t+Lag}, soi(t)} for positive values of Lag ... or ... shouldn't the title be soi & rec?? oops.

Now, try this

ccf(soi,rec)

and you get

ccf2

What you're seeing is corr{soi(t+Lag), rec(t)} versus Lag. So on the positive side of Lag, rec leads, and on the negative side of Lag, soi leads.

We're not done with this yet. If you want to do a regression of x on lag(x,-1), for example, you have to "tie" them together first. Here's an example.

x = arima.sim(list(order=c(1,0,0), ar=.9), n=20)  # 20 obs from an AR(1)
xb1 = lag(x,-1)
## you wouldn't regress x on lag(x,1) because that would be progress :)

##-- WRONG:
cor(x,xb1) # correlation
[1] 1 ... is one?
lm(x~xb1) # regression
Coefficients:
(Intercept) xb1
6.112e-17 1.000e+00

do it the WRONG way and you'll think x(t)=x(t-1)


##-- RIGHT:
y=cbind(x,xb1) # tie them together first
lm(y[,1]~y[,2]) # regression
Coefficients:
(Intercept) y[, 2]
0.5022 0.7315


##-- OR:
y=ts.intersect(x,xb1) # tie them together first this way
lm(y[,1]~y[,2]) # regression
Coefficients:
(Intercept) y[, 2]
0.5022 0.7315

cor(y[,1],y[,2]) # correlation
[1] 0.842086


By the way, (Intercept) is used correctly here.

R does warn you about this (type ?lm and scroll down to "Using time series"), so consider this a heads-up, rather than an issue. See our little tutorial for more info on this.

ISSUE 6: Why do you want to see the zero lag value of the ACF, and what's this mess?

When you're trying to fit an ARMA model to data, one of the first things you do is look at the ACF and PACF of the data. Let's try this for a simulated MA(1) process. Here's how:

MA1=arima.sim(list(order=c(0,0,1), ma=.5), n=100)
par(mfcol=c(2,1))
acf(MA1,20)
pacf(MA1,20)
and here's the output:

ma1

What's wrong with this picture? First, the two graphs are on different scales. The ACF axis goes from -.2 to 1, whereas the PACF axis goes from -.2 to .4. Also, the lag axis on the ACF plot starts at 0 (the 0 lag ACF is always 1 so you have to ignore it or put your thumb over it), whereas the lag axis on the PACF plot starts at 1.

So, instead of getting a nice picture by default, you get a messy picture. Ok, the remedy is as follows:

par(mfrow=c(2,1))
acf(MA1,20,xlim=c(1,20)) # set the x-axis limits to start at 1 then
# look at the graph and note the y-axis limits
pacf(MA1,20,ylim=c(-.2,1)) # then use those limits here
and here's the output:

ma1a

Looks nice, but who wants to get carpal tunnel syndrome sooner than necessary? Not me. So I wrote an R function called acf2 that will do everything at once and save you some time and typing. You can get acf2 on the web page for the text under R CODE (Ch 1-5) - use the blue bar on top of this page to get there.

0
Importantes notas para tener en cuenta. (Stoffer)

ISSUE 1: When is the intercept the mean?

When fitting ARIMA models, R calls the estimate of the mean, the estimate of the intercept. This is ok if there's no AR term, but not if there is an AR term. For example, suppose x(t) = α + φ*x(t-1) + w(t) is stationary. Then taking expectations we have μ = α + φ*μ or α = μ*(1-φ). So, the intercept, α, is not the mean, μ, unless φ=0. In general, the mean and the intercept are the same only when there is no AR term. Here's a numerical example:

# generate an AR(1) with mean 50
set.seed(66) # so you can reproduce these results
x = arima.sim(list(order=c(1,0,0), ar=.9), n=100) + 50
mean(x)
[1] 50.60668 # the sample mean is close

arima(x, order = c(1, 0, 0))
Coefficients:
ar1 intercept <-- here's the problem
0.8971 50.6304 <-- or here, one of these has to change
s.e. 0.0409 0.8365

The result is telling you the estimated model is
x(t) = 50.6304 + .8971*x(t-1) + w(t)
whereas, it should be telling you the estimated model is
x(t)−50.6304 = .8971*[x(t-1)−50.6304] + w(t)
or
x(t) = 5.21 + .8971*x(t-1) + w(t). Note that 5.21 = 50.6304*(1-.8971).

The easy thing to do is simply change "intercept" to "mean":

  Coefficients:
ar1 mean
0.8971 50.6304
s.e. 0.0409 0.8365


I should mention that I reported this flub to the R folks, but I was told that it is a matter of opinion. But, if you use ar.ols, you get anotheR opinion:
ar.ols(x, order=1, demean=F, intercept=T)
Coefficients:
1
0.9052 <-- estimate of φ
Intercept: 4.806 (2.167) <-- yes, it IS the intercept as you know and love it

Note that arima() uses MLE, whereas ar.ols() uses OLS to fit the model, and hence the differences in the estimates. One thing is certain, the use of the term intercept in R is open to interpretation, which is not exactly an optimal situation.

ISSUE 2: Why does arima fit different models for different orders? And what does xregdo?

When fitting ARIMA models with R, a constant term is NOT included in the model if there is any differencing. The best R will do by default is fit a mean if there is no differencing [type ?arima for details]. What's wrong with this? Well (with a time series in x), for example:

arima(x, order = c(1, 1, 0))          #(1)

will not produce the same result as

arima(diff(x), order = c(1, 0, 0))    #(2)

because in (1), R will fit the model [with ∇x(s) = x(s)-x(s-1)]
∇x(t)= φ*∇x(t-1) + w(t)
whereas in (2), R will fit the model
∇x(t) = α + φ*∇x(t-1) + w(t).

If there's drift (i.e., α is NOT zero), the two fits can be extremely different and using (1) will lead to an incorrect fit and consequently bad forecasts (see Issue 3 below).

If α is NOT zero, then what you have to do to correct (1) is use xreg as follows:

arima(x, order = c(1, 1, 0), xreg=1:length(x))    #(1+)

Why does this work? In symbols, xreg = t and consequently, R will replace x(t) with x(t)-β*t; that is, it will fit the model
∇[x(t) - β*t] = φ*∇[x(t-1) - β*(t-1)] + w(t).
Simplifying,
∇x(t) = α + φ*∇x(t-1) + w(t) where α = β*(1-φ).

If you want to see the differences, generate a random walk with drift and try to fit an ARIMA(1,1,0) model to it. Here's how:

set.seed(1)           # so you can reproduce the results
v = rnorm(100,1,1) # v contains 100 iid N(1,1) variates
x = cumsum(v) # x is a random walk with drift = 1
plot.ts(x) # pretty picture...

arima(x, order = c(1, 1, 0)) #(1)

Coefficients:
ar1
0.6031
s.e. 0.0793

arima(diff(x), order = c(1, 0, 0)) #(2)

Coefficients:
ar1 intercept <-- remember, this is the mean of diff(x)
-0.0031 1.1163 and NOT the intercept
s.e. 0.1002 0.0897

arima(x, order = c(1, 1, 0), xreg=1:length(x)) #(1+)

Coefficients:
ar1 1:length(x) <-- this is the intercept of the model
-0.0031 1.1169 for diff(x)... got a headache?
s.e. 0.1002 0.0897

Let me explain what's going on here. The model generating the data is
x(t) = 1 + x(t-1) + w(t)
where w(t) is N(0,1) noise. Another way to write this is
[x(t)-x(t-1)] = 1 + 0*[x(t-1)-x(t-2)] + w(t)
or
∇x(t) = 1 + 0*∇x(t-1) + w(t)
so, if you fit an AR(1) to ∇x(t), the estimates should be, approximately, ar1 = 0 and intercept = 1.

Note that (1) gives the WRONG answer because it's forcing the regression to go through the origin. But, (2) and (1+) give the correct answers expressed in two different ways.

ISSUE 3: Why does predict.Arima give strange forecasts?

If you want to get predictions from an ARIMA(p,d,q) fit when there is differencing (i.e., d > 0), then the previous issue continues to be a problem. Here's an example using the global temperature data from Chapter 3. In what you'll see below, the first method gives the wrong results and the second method gives the correct results. I think it's just best to stay away from the first method. If you use sarima and sarima.for, then you'll avoid these problems.

 u=read.table("/mydata/globtemp2.dat")  # read the data
gtemp=ts(u[,2],start=1880,freq=1) # yearly temp in col 2
fit1=arima(gtemp, order=c(1,1,1))
fore1=predict(fit1, 15)
nobs=length(gtemp)
fit2=arima(gtemp, order=c(1,1,1), xreg=1:nobs)
fore2=predict(fit2, 15, newxreg=(nobs+1):(nobs+15))
par(mfrow=c(2,1))
ts.plot(gtemp,fore1$pred,col=1:2,main="WRONG")
ts.plot(gtemp,fore2$pred,col=1:2,main="RIGHT")

Here's the graphic:

forecasting

ISSUE 4: tsdiag.Arima gives the wrong p-values

If you use tsdiag() for diagnostics after an ARIMA fit, you will get a graphic that looks like this:

tsdiag
The p-values shown for the Ljung-Box statistic plot are incorrect because the degrees of freedom used to calculate the p-values are lag instead of lag - (p+q). That is, the procedure being used does NOT take into account the fact that the residuals are from a fitted model. And YES, at least one R core developer knows this.

0
Granger: casualidad, causalidad y cointegración - La Nueva España - Diario Independiente de Asturias

Granger: casualidad, causalidad y cointegración - La Nueva España - Diario Independiente de Asturias: "- Enviado mediante la barra Google"

La clave está, pues, en ver si hay o no hay «cointegración», un término que se inventó en 1981 para referirse a aquellas series que siendo «malas» (no estacionarias) terminan siendo buenas (estacionarias) cuando se juntan, se combinan. Aunque no es habitual que cuando se juntan los malos el resultado sea bueno, lo cierto es que la vida nos da sorpresas y a veces ocurre.

0
Cointegración

IN: Cointegración.
La econometría de series temporales se encuentra con un problema al medir las
relaciones entre aquellas variables que tienen una tendencia temporal. Este problema
puede llegar a que se consideren significativas relaciones completamente espurias.
Las variables que tienen una tendencia temporal definida se denominan “no
estacionarias”. Las estimaciones de regresiones con variables no estacionarias son
espurias salvo que estas estén cointegradas. Dos variables no estacionarias cointegradas
son aquellas cuyos residuos son estacionarios. Si los residuos son estacionarios las
estimaciones de variables no estacionarias son superconsistentes.

IN: Cointegración: "Econometría" . ¿qué es la cointegración?
Muy sencillo: supongamos que dos series temporales, xt e yt, son estacionarias de orden 1 (es decir, son I(1); para órdenes superiores de integración como I(2), I(3), etc. el problema se complicaría algo más). Se dice que dichas variables están cointegradas cuando puede practicarse una regresión lineal del siguiente tipo:

yt = a·xt + ut


De tal forma que los residuos (errores de ajuste) de la regresión, ut = yt – a·xt sean I(0), esto es, estacionarios. Por tanto, en su versión más sencilla, la cointegración exige que se verifiquen dos condiciones básicas:

  • Que dos variables sean integradas de orden 1.
  • Que exista una combinación lineal de ambas que sea estacionaria de orden 0.

El concepto de cointegración es relativamente reciente. Fue acuñado en 1987 por C.W.J. Granger, reconocido economista británico que falleció el año pasado y que fue premio Nobel de Economía en 2003 junto con su inseparable compañero Robert Engle (sí, el mismo de los modelos ARCH). En sus investigaciones, Granger observó que la mayoría de los economistas utilizaban series no integradas para estimar relaciones entre ellas, lo cual podía conducir a obtener relaciones espurias, es decir, que se diera el caso de que dos variables estuvieran aparentemente relacionadas cuando en realidad no lo estaban, existiendo una tercera variable desconocida que las relacionaba. Tal sería el caso por ejemplo de un investigador que detecta una relación entre estatura e inteligencia en un colegio. No es que por ser más alto un niño sea más inteligente, sino que al ser más maduros los niños de mayor edad (y, por tanto, mayor estatura) las habilidades y destrezas cognitivas también son mayores. Pues bien, cuando las variables estaban cointegradas este problema desaparece ya que los residuos obtenidos son I(0) indicando que no queda ninguna estructura pendiente de modelizar.

Paquete URCA para realizar tests de cointegración
http://cran.r-project.org/web/packages/urca/index.html
http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=urca:plot-methods

cointegracion

0
Estacionariedad y Cointegración. Teoría. Aplicación sobre series climáticas.

Notas Sobre Analisis de Series de Tiempo


Aplicación sobre series climáticas

Temp Analysis2

lunes, 10 de enero de 2011

0
Time Series Regression of Global Temperature, El Nino – LaNina, and Volcanic Events | Climate Charts & Graphs

0
ARIMA. Modelos no estacionarios estacionales


Muchas series temporales no son estacionarias, ya sea porque presentan tendencias o por efectos estacionales. En particular, las caminatas aleatorias, que caracterizan muchos tipos de series, no son estacionarias pero pueden transformarse en series estacionarias utilizando la diferenciación de primer orden. Como la serie diferenciada necesita ser agregada (o integrada) para recuperar la serie original, el proceso estocástico subyacente se llama media móvil integrada autorregersiva (ARIMA).
Los procesos ARIMA se pueden extender para incluir términos estacionales, dando paso a modelos ARIMA estacionales no-estacionarios (SARIMA).
Las series no-estacionarias porque la varianza está correlacionada en serie (heterocedástica condicional) resultan en períodos de volatilidad (un cambio claro en la varianza). Una aproximación para modelar este tipo de series es el modelo heterocedástico condicional autorregresivo (ARCH). Asimismo, podemos utilizar la generalización del modelo ARCH (GARCH).
  • Modelos ARIMA no-estacionales
Una serie temporal {X(t)} sigue un proceso ARIMA(p,d,q) si la diferencia dth de la serie es un proceso ARMA(p,q). Es decir, sea Y(t)=(1-B)^d*X(t) entonces theta_p(B)*Y(t)=phi_q(B)*W(t). Por tanto, podemos escribir el proceso ARIMA como_ thet_p(B)(1-B)^d*X(t)=phi_q(B)*W(t), donde theta_p y phi_q son polinomios de orden p y q, respectivamente..

Ejemplos:
  • ARIMA(0,1,1) o IMA(1,1): X(t)=X(t-1)+W(t)+beta*W(t-1), donde beta es el parámetro del modelo, o (1-B)*X(t)=(1+beta*B)*W(t).
  • ARIMA(1,1,0) o ARI(1,1): X(t)=alfa*X(t-1)+X(t-1)-alfa*X(t-2)+W(t) donde alfa es el parámetro modelado, o (1-alfa*B)(1-B)*X(t)=W(t).

Aplicación en R: ajuste ARIMA(1,1,1) para series simuladas
    #Ajuste del ARIMA(1,1,1) para series simuladas mediante el modelo xt = 0.5xt-1+xt-1-0.5xt-2+wt+0.3wt-1.
#manualmente
set.seed(1)
x <- w <- rnorm(1000); for (i in 3:1000) x[i] <- 0.5 * x[i - 1] + x[i - 1] - 0.5 *x[i - 2] + w[i] + 0.3 * w[i - 1]
arima(x, order = c(1, 1, 1))
#o automáticamente
x <- arima.sim(model = list(order = c(1, 1, 1), ar = 0.5,ma = 0.3), n = 1000)
arima(x, order = c(1, 1, 1))



Aplicación en R: ajuste ARIMA(0,1,1) para la producción de cerveza australiana
    #Ajuste del IMA(1,1) o ARIMA(0,1,1) para la producción de cerveza australiana
www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, he = T)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Beer.ima <- arima(Beer.ts, order = c(0, 1, 1)) #la serie tiene una tendencia creciente, por loque ajustamos un modelo IMA(1,1) (representa una tendencia linal con un ruido blanco)
Beer.ima #parámetros estimados xt=Xt-1+Wt-.33*wt-1.
acf(resid(Beer.ima)) #corelograma de los residuales. tiene picos en cada año y sugiere que se necesita un término estacional.
Beer.1991 <- predict(Beer.ima, n.ahead = 12) #predecimos los valores para los próximos 12 años
sum(Beer.1991$pred) #producción anual total para el año 1991



Aplicación en R: diferenciación de la serie de electricidad
#Ajuste de ARIMA no estacional a la serie de electricidad. Diferenciación para eliminar tendencia.
www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, he = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
layout(c(1, 1, 2, 3))
plot(Elec.ts)
plot(diff(Elec.ts))
plot(diff(log(Elec.ts)) #eliminamos la tendencia creciente




  • Modelos ARIMA estacionales
El modelo ARIMA estacional usa la diferenciación en lag igual al número de estaciones (s) para elimina el efecto estacional aditivo. Al igual que la diferenciación a lag=1 que remueve la tendencia, la diferenciación a lag=s introduce un término de medias móviles, además del término autorregresivo en lag=s. El ARIMA estacional ARIMA(p,d,q)(P,D,Q)s puede expresarse usando el operador de retardo B:

Theta_P(B^s)*theta_p(B)(1-B^s)^D(1-B)^d*x(t)=Phi_Q(B^s)*phi_q(B)*W(t)

donde Theta_P, theta_p, Phi_Q y phi_q son polinomios de orden P, p, Q y q, respectivamente. El modelo es no-estacionario, pero si D=d=0 y las raíces de la ecuación característica son mayores a la unidad, el modelo resultante será estacionario.

Ejemplos:
  • AR(1) con período estacional de 12 unidades o ARIMA(0,0,0)(1,0,0)12: es el modelo xt = alfa*xt−12 + wt. Este modelo es apropiado para datos mensuales donde mes actual solo está influenciado por el valor del mes del año anterior.
  • ARIMA(0,1,0)(1,0,0)12 o xt = xt−1 + alfa*xt−12 – alfa*xt−13 + wt, es un modelo con tendencia estocástica e influencias estacionales. El cambio en el tiempo t depende del cambio en el mismo tiempo (i.e. mes) del año anterior.
  • MA(1) con período estacional de 4 unidades o ARIMA(0,0,0)(0,0,1)4: es el modelo xt=wt-beta*wt-4, estacionario y solo apto para datos sin tendencia.
  • ARIMA(0,1,0)(0,0,1)4 o xt=xt-1+wt-beta*wt-4: para datos con tendencia estocástica e influencias estacionales cuatrimestrales.
  • ARIMA(0,0,0)(0,1,1)4 o xt = xt−4 + wt – beta*wt−4: para datos donde el término estacional contiene una tendencia estocástica (la diferenciación se aplica al período estacional). con término el modelo se puede extender

Criterios de selección de modelos ARIMA:
Los modelos SARIMA pueden tener un gran número de parámetros y combinaciones de ellos. Por tanto, será apropiado considerar un amplio rango de modelos posibles y elegir entre ellos según un criterio apropiado (e.g. AIC). Para ello deberemos tener cuidado con el sobre-ajuste (el cual produce mejores AIC). Una vez que seleccionamos el mejor modelo a ajustar, realizamos el correlograma de los residuos para verificar que sigue un comportamiento de ruido blanco.

Aplicación en R: selección entre ajuste ARIMA(1,1,0)(1,0,0)12 y ARIMA(0,1,1)(0,0,1)12 para series de producción eléctrica, luego buscamos el mejor modelo según un procedimiento automático y lo analizamos.

    #Ajuste de ARIMA estacional a la serie de electricidad.
#ajustamos dos modelos ajuste ARIMA(1,1,0)(1,0,0)12 y ARIMA(0,1,1)(0,0,1)12 y evaluamos cuál se ajusta mejor a los datos mediante el criterio AIC
AIC (arima(log(Elec.ts), order = c(1,1,0),seas = list(order = c(1,0,0), 12)))
AIC (arima(log(Elec.ts), order = c(0,1,1),seas = list(order = c(0,0,1), 12)))
#el modelo ARI estacional provee el mejor ajuste, con el menor AIC
#para chequear varios modelos podemos utilizar la siguiente función con el criterio CAIC (que permite eliminar la sobre-parametrización)
get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1)) #construimos una función para obtener el mejor modelo arima según el criterio CAIC
{
best.aic <- 1e8
n <- length(x.ts)
for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
{
fit <- arima(x.ts, order = c(p,d,q),
seas = list(order = c(P,D,Q),
frequency(x.ts)), method = "CSS")
fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
if (fit.aic < best.aic)
{
best.aic <- fit.aic
best.fit <- fit
best.model <- c(p,d,q,P,D,Q)
}
}
list(best.aic, best.fit, best.model)
}
best.arima.elec <- get.best.arima( log(Elec.ts),maxord = c(2,2,2,2,2,2)) #observamos que el mejor modelo ajustado usando términos de orden mayor a 2, es ARIMA(0,1,1)(2,0,2)12.
best.fit.elec <- best.arima.elec[[2]]
acf( resid(best.fit.elec) ) #los residuales parecen aproximarse a un ruido blanco
best.arima.elec [[3]]
ts.plot( cbind( window(Elec.ts,start = 1981),exp(predict(best.fit.elec,12)$pred) ), lty = 1:2) #se utiliza un factor de correción sesgado para los valores predichos. Sin embargo, esa corrección parece innecesaria en este caso porque la desviación estandar de los residuos es pequeña en comparación con las predicciones.



RESUMEN DE PROCESOS ARIMA ESTACIONALES (SARIMA)

  • Procesos no estacionarios, donde analizaremos la falta de estacionariedad en la media (comportamiento estacional).
  • Generalmente podemos incorporar la estacionalidad dentro de un modelo ARIMA de una forma multiplicativa, resultando en un modelo ARIMA estacional multiplicativo.
  • Concepto de estacionalidad y sus tipos:
  1. Una serie estacional presenta valores que no son constantes pero varían con una pauta cíclica. Cuando E[Xt]=E[Xt+s] decimos que s es el período de estacionalidad, y éste define el número de observaciones que forman el ciclo estacional: s=12 (serie mensual), s=4 (trimestral), s=7 (semanal).
  2. El modelo más simple trata un efecto constante que se suma a los valores de la serie: Xt=St+nt. La serie se escribe como suma de un componente estacional St y un proceso estacionario nt (con mu su media).
  3. El componente estacional St puede tener comportamiento: 1) determinista (función constante para el mismo mes en distintos años) St=St+ks, 2) estacionario (evoluciona en el tiempo y su evolución es estacionaria, oscilando alrededor de un valor medio) St=mu+vt, donde vt es un proceso estacionario de media cero que introduce variabilidad en el año, o 3= no-estacionario (es cambiante sin ningún valor medio fijo) St=St-s+vt.
En los 3 casos podemos convertir una serie estacional en estacionaria aplicando una diferencia estacional. Definimos el operador diferencia estacional de período s: diffs=1-B^s, es decir diffsXt=(1-B^s)Xt=Xt-Xt-s p diffsXt=diffsSt+diffsnt para nuestro modelo más simple de estacionalidad.
En el modelo ARIMA estacional podemos convertir series no estacionarias en estacionarias tomando diferencias regulares (entre períodos consecutivos) o eliminar la estacionalidad mediante diferencias estacionales, ambos mediante la transformación: Wt=diffsDdiffdXt, donde D es el número de diferencias estacionales (generalmente D=1) y del número de diferencias regulares (d<=3). Cuando existe dependencia estacional podemos generalizar el modelo ARMA para series estacionarias incorporando además la dependencia regular (asociada a intervalos de medida de la serie) y la dependencia estacional (asociada a observaciones separadas por s períodos). Debemos modelar 2 tipos de dependencias:
  1. Incorporar la dependencia estacional a la regular, añadiendo a los operadores AR o MA el operador B, términos B^s para representar la dependencia entre observaciones separadas por s períodos.
  2. Modelar de forma separada la dependencia regular y la estacional, y construir el modelo incorporando ambas de forma multiplicativa: modelo ARIMA estacional multiplicativo. Esta opción es la más simple.
Esta clase de modelos se escriben de forma simplificada como el modelo ARIMA(P,D,Q)sx(p,d,q).
Nota: el modelo ARIMA estacional multiplicativo se basa en la hipótesis central de que la relación de dependencia estacional (modelo estacional) es la misma para todos los períodos.

En el ACF:
fijarse en los retardos iniciales 1,2,3,4,… para identificar la estructura regular, y en los retardos estacionales s, 2s, 3s,… para identificar la estructura estacional. La interacción alrededor de los coeficientes estacionales puede entonces utilizarse como confirmación de la identificación realizada.
En el PACF: fijarse en los retardos iniciales 1,2,3,4… para identificar la estructura regular y en los retardos estacionales s,Ds,3s,… para la identificación de la estructura estacional.


IDENTIFICACIÓN DE LOS POSIBLES MODELOS ARIMA

  1. Selección del modelo de la serie estacionaria: identificar la estructura no-estacionaria (si existe) y después la estructura ARMA estacionaria.
  • a. Decidir qué transformaciones aplicar para convertir la serie observada en una serie estacionaria. Transformar la serie para que tenga varianza constante y/o diferenciar la serie para que tenga media constante.
  • -Si la impresión visual de la serie es que la variabilidad no aumenta con el tiempo sino con el nivel de la serie, nos fijaremos en el gráfico de desviación típica (medida de la variabilidad) vs. Media local (medida del nivel) para confirmar o no nuestra sospecha visual. Cuando la variabilidad aumenta linealmente con el nivel de la serie, tomaremos logaritmos para que la variabilidad sea constante. En general, aplicamos la transformación Box-Cox Yt={xt^(1-alfa)-1}/[1-alfa], donde alfa lo estimamos como la pendiente del gráfico log(si) vs. Log(mean(xi)), donde si es la desviación típica y mean(xi) es la media. Excepciones: si la varianza cambia sin relación con el nivel o si es constante.
  • b. Determinar un modelo ARMA para la serie estacionaria (p y q de su estructura AR y MA) y si el proceso es estacional, determinar los órdenes P, Q de la estructura ARMA estacional.
  • Estimaremos todos los modelos que consideremos posibles para explicar la serie y después elegir entre ellos.

  1. Estimación: los parámetros AR y MA del modelo se estiman por máxima verosimilitud y se obtienen sus errores estándar y los residuos del modelo.
  2. Diagnosis: se comprueba que los residuos no tienen estructura de dependencia y siguen un proceso de ruido blanco. Si los residuos tienen estructura modificaremos el modelo para incorporarla y repetiremos las 3 etapas anteriores hasta obtener un modelo adecuado.

Referencias
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). Introductory Time Series with R. Springer-Verlag.

0
ARMA. Modelos estacionarios no estacionales

Los modelos de regresión pueden dar cuenta de componentes no estacionarios, sin embargo, no consideran la correlación temporal de los residuales. Por ello, consideraremos modelos estacionarios que pueden ser utilizados para series de residuos que contengan tendencias o ciclos estacionales no obvios. El modelo estacionario ajustado puede combinarse con un modelo de regresión para mejorar los pronósticos.


1. Modelos de media móvil MA(q)

Definición: un proceso MA de orden 1 es una combinación lineal del término de ruido blanco actual y los q términos de ruido blanco del pasado más reciente.

X(t) = w(t) +beta_1*w(t−1) + . . . + beta_q*w(t−q)

donde {w(t)} es un ruido blanco con media cero y varianza sigma^2.
También podemos escribir esta ecuación en términos del operador de cambio de retardo B:

X(t)=(1+beta_1*B+…+beta_q*B^q)*w(t)=phi_q(B)*w(t)

donde phi_q es un polinomio de orden q.
Como el proceso MA consiste en la suma de términos de ruido blanco estacionarios, son estacionarios y tienen media y autocovarianza invariantes en el tiempo.

Función “arima”: El modelo MA(q) puede ser ajustado a los datos utilizando la función “arima” seleccionando el orden de los parámetros como c(0,0,q). Por defecto, la función “arima” no resta la media y estima un término de intercepto. Asimismo, minimiza la suma condicional de cuadrados para estimar los valores de los parámetros (para verlos debemos seleccionar “method”=c(“CSS”) o utilizarlos como entradas para estimación por máxima verosimilitud). Podemos seleccionar el valor de la media como cero, en logar de estimar el intercepto, utilizando “incluse.mean=FALSE” dentro de la función “arima”.

Aplicación en R: ajuste de modelos MA a series simuladas según el modelo MA(3)
set.seed(1)

b <- c(0.8, 0.6, 0.4)
x <- w <- rnorm(1000)
for (t in 4:1000) {
for (j in 1:3)
{
x[t] <- x[t] + b[j] * w[t - j]
}
plot(x, type = "l")
acf(x)
x.ma <- arima(x, order = c(0, 0, 3))
x.ma #observamos los parámetros estimados y los IC al 95%, que contienen los valores de los parámetros subyacentes.



Aplicación en R: ajuste del modelo MA(1) para la serie de tasas de intercambio
#Ajuste de modelo MA(1). Serie de tasa de intercambio.

www <- "http://www.massey.ac.nz/~pscowper/ts/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1)) #ajuste MA(1)
x.ma
acf(x.ma$res[-1]) #el modelo no provee un ajuste satisfactorio, los residuos no son una realización realista del ruido blanco



2. Modelos mixtos ARMA


Definición: cuando juntamos los términos AR y MA en una única expresión obtenemos un modelo de media móvil autorregresivo (ARMA) de orden (p,q):

X(t)=alfa_1*X(t-1)+…+ alfa_p*X(t-p)+W(t)+beta_1*W(t-1)+…+beta_q*W(t-q)

donde {W(t)} es un ruido blanco. En términos del operador de retardos B, tenemos:

theta_p(B)*X(t)=phi_q(B)*W(t)

El modelo ARMA(p,q) tiene las siguientes características:
• Es estacionario cuando las raíces de theta son mayores a la unidad, en valor absoluto.
• Es invertible cuando las raíces de phi son mayores a la unidad, en valor absoluto.
• El modelo AR(p) es un caso especial del ARMA(p,0)
• El modelo MA(q) es un caso especial del ARMA(0,q)
Parsimonia sobre los parámetros: el modelo ARMA es más eficiente (i.e. requiere menores parámetros) que los modelos AR o MA solos.
Redundancia de los parámetros: cuando theta y phi tienen un factor común, el modelo estacionario puede simplificarse.

Función “arima.sim”: los procesos ARMA y ARIMA pueden simularse con esta función.

Aplicación en R: ajuste una series simulada según el modelo ARMA(1,1) con alfa=-.6 y beta=.5.
   #Ajuste de modelo ARMA(1,1) para una serie simulada con el modelo ARMA(1,1) y parámetros alfa=-.6 t beta=.5.

set.seed(1)
x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1))) #como era de esperar las estimaciones muestrales para alfa y beta son cercanos a los parámetros del modelo subyacente


Aplicación en R: ajuste de modelos ARMA para las series de tasas de intercambio.
#Ajuste de modelos ARMA para las series de tasa de intercambio.

x.ma <- arima(x.ts, order = c(0, 0, 1)) #ajuste MA(1) o ARIMA(0,0,1)
x.ar <- arima(x.ts, order = c(1, 0, 0)) #ajuste AR(1) o ARIMA(1,0,0)
x.arma <- arima(x.ts, order = c(1, 0, 1)) #ajuste ARMA(1,1) o ARIMA(1,0,1)
AIC(x.ma); AIC(x.ar) ; AIC(x.arma) #comparación de los modelos utilizando el criterio AIC. El ARMA(1,1) da mejor ajuste a los datos, seguido por el AR(1) y finalmente el MA(1).
x.arma
acf(resid(x.arma)) #los residuales del ARMA(1,1) tienen pequeñas autocorrelaciones, lo cual es consistente con la realización de un ruido blanco.


Aplicación en R: Ajuste de modelos ARMA para las series de producción eléctrica.
#Ajuste de modelo ARMA para la producción de electricidad.

www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
tsdisplay(Elec.ts) #existe una tendencia positiva y un ciclo estacional regular. La varianza aumenta con el tiempo sugiriendo una transformación log.
Time <- 1:length(Elec.ts)
Imth <- cycle(Elec.ts)
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth)) #se ajuste un modelo de regersión a los logaritmos de la serie original.
acf(resid(Elec.lm)) #el correlograma de los residuales parece ciclar con período de 12 meses sugiriendo que las variables indicadoras mensuales no son suficientes para dar cuenta de la estacionalidad en la serie.
#podemos dar cuenta de esto mediante un modelo no-estacionario con un componente estacional estocástico
best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:2)
{
for (j in 0:2)
{
fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0,j)))
if (fit.aic < best.aic)
{
best.order <- c(i, 0, j)
best.arma <- arima(resid(Elec.lm), order = best.order)
best.aic <- fit.aic
}
}
}
best.order
acf(resid(best.arma)) #correlograma del modelo ARMA de los residuales con mayor AIC

#predecimos la producción eléctrica para cada mes de los siguientes 3 años
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12,3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991,freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)


Aplicación en R: ajuste de modelos ARMA para los datos de ondas en experimentos de tanques/baldes.
 #Ajuste ARMA. Datos del tanque de ondas.

www <- "http://www.massey.ac.nz/~pscowper/ts/wave.dat"
wave.dat <- read.table(www, header = T)
attach (wave.dat)
layout(1:3)
plot (as.ts(waveht), ylab = 'Wave height')
acf (waveht); pacf (waveht) #sugiere que p debe ser al menos 2
wave.arma <- arima(waveht, order = c(4,0,4)) #el mejor modelo aRMA(p,q) obtenido según la varianza mínima de los residuales, se obtiene para p=q=4. Ajustamos un ARMA(4,4)
acf (wave.arma$res[-(1:4)]); pacf (wave.arma$res[-(1:4)]) #se adecua a un proceso de ruido blanco
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')



RESUMEN DE PROCESOS ARMA
  • Son combinación de las propiedades de los procesos AR y MA y permiten representar de forma “escueta” (utilizan pocos parámetros) procesos cuyas primeros q coeficientes son cualesquiera, mientras que los siguientes decrecerán según leyes simples.
ARMA(1,1)
Xt=phi*Xt-1+epsilont-phi1*epsilont-1 o (1-phi1B)Xt=(1-theta1B)epsilon, donde |phi1|<1>
ARMA(p,q)

(1-phi1B-…-phipB^p)Xt=(1-theta1B-…-thetaqB^q)epsilont o phip(B)Xt=thetaq(B)epsilont.



Referencias
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). Introductory Time Series with R. Springer-Verlag.
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). http://elena.aut.ac.nz/~pcowpert/ts/