Construyendo un modelo ARMA para series de tiempo

Sobre los errores o residuos

En las entradas previas no mencioné mucho respecto a los residuos de los ajustes o modelos. Lo que hice fue mostrar algunas de las gráficas más usuales y mostrar como se pueden obtener en R project.

Antes de mencionar algo sobre los errores, hago un ejemplo sencillo de lo que es un ruido blanco. No menciono la definición formal pero trato de explicar cual es la idea.

#Ejemplo Ruido Blanco
#Generamos los datos

ruido_blanco=rnorm(1000,0,1)

#Graficamos la serie de tiempo

plot.ts(ruido_blanco, main="Ejemplo de Ruido Blanco", xlab="Tiempo", ylab="Valores",col="6")

#Graficamos su correlograma

acf(ruido_blanco,main="Correlograma",col="2",lag=50)

#Calculamos la media

mean(ruido_blanco)
#[1] 0.03784017

Ejemplo_RB

Correlograma_RB

Un ruido blanco es una serie tal que su media es cero, la varianza es constante y es incorrelacionada. Lo que se hice fue simular valores de una serie de tiempo proveniente de una distribución normal con media cero y varianza uno. Se calculé la media con el comando mean y se analizo la incorrelación por medio del correlograma generado de la función afc . Se observa en la gráfica que sus valores para el t>0 son casi cero o que están debajo de las bandas azules de la gráfica (bandas de Bartlett) .

Hay dos cosas importantes en la definición de ruido blanco, pedir que la varianza sea constante con respecto al tiempo es algo bastante relevante. Esa propiedad se llama homocedasticidad. Cuando se analiza la regresión lineal se pueden hace varias pruebas estadísticas al respecto, pero con las series de tiempo basta analizar el comportamiento de la función de autocorrelación que se gráfica en los correlogramas.

También se puede probar que una serie es un ruido blanco por medio de prueba de hipótesis. Para esto se puede hacer uso de  dos posibles pruebas Ljung-Box o Durbin-Watson.

#Ejemplo Ruido Blanco
#Prueba Ljung-Box

Box.test(ruido_blanco)

La relevancia de los ruidos blancos es que son piezas fundamentales de la construcción de los modelos arma. Por eso es importante tener ubicadas sus propiedades y como simularlos.

Cabe mencionar que existe una relación entre las caminatas aleatorias y los ruidos blancos. En la entrada anterior, ¿qué es una serie de tiempo?, puse dos ejemplos de caminatas aleatorias y ruidos blancos en una misma gráfica. En resumen, se dice que la diferenciación de una caminata aleatoria es un ruido blanco. Explico más adelante a que se refiere esto.

Algo de procesos estacionarios

No doy la definición formal de procesos estocástico estacionario, pero la razón de este concepto con las series de tiempo tienen que ver con le hecho de pensar la serie de tiempo Y(t) como una familia de variables aleatorias que tienen una distribución conjunta y que el valor que observamos de la seria es el valor particular de las variables.

La idea es que la serie es un caso particular de valores posibles de las variables aleatorias y por lo cual se puede definir p(Y1,Y2,…,Yt) y p(Y(t+1)|Y1,Y2,…,Yt). Es decir, su probabilidad conjunta y su probabilidad condicional. Ser estacionario tiene que ver con la posibilidad de cambiar t y no afectar la distribución conjunta ni la distribución condicional.

Quizás es enredoso pensar en el concepto de que la seria sea estacionaria, pero es la idea teórica que permite desarrollar las técnicas de análisis.

En resumen, lo que termina significando el concepto es que deseamos que la media de la serie sea constante, que la varianza también y que la covarianza solo dependa de la distancia temporal entre los datos. Estos supuesto son parecidos a los que mencioné de los ruidos blancos.

Simplificando el concepto se puede definir una serie estacionaria en covarianza, con lo cual solo se pide que la media sea constante y que la covarianza sea  solo dependiente de la distancia temporal entre las variables.

Todo lo teórico antes mencionado se reduce analizar la función de autocorrelación o el correlograma. Si la gráfica del correlograma muestra que las barras tienden a cero o va teniendo valores muy pequeños, podemos considerar que la serie que analizamos es estacionaria por lo menos en covarianza.

#Ejemplo de series estacionarias y no estacionarias
#Simulados un ruido blanco gaussiano

ruido_blancog=rnorm(1000,0,1)
acf(ruido_blancog,lag=200,col="2",main="Ejemplo de serie estacionaria",xlab="Tiempo",ylab="Valores")

#Tomamos datos de las tasa de interés de USA
rate=read.table("m-gs10.txt",head=T)

rate=ts(rate[,4])

acf(rate,lag=200,col="2",main="Ejemplo de serie no estacionaria",xlab="Tiempo",ylab="Valores")

Ejemplo_ts_estacionaria

Ejemplo_ts_no-estacionaria

Muchas series en la práctica no son estacionarias ni en covarianza. Lo que se hace en esos casos  es calcular la diferencias la serie o transformarla por alguna función y luego diferenciar. El proceso anterior se define como proceso homogéneo no estacionario y el grado de diferenciación es el orden de homogeneidad. Esto lo muestro con unos ejemplos en el siguiente código.

#Ejemplos de diferenciar la serie y de transformar y diferencias

rate=read.table("m-gs10.txt",head=T)
rate=ts(rate[,4])

plot(diff(rate,lag=1),main="Ejemplo de Diferenciación", xlab="Tiempo", ylab="Valores", col="2")

acf(diff(rate,lag=1),main="Ejemplo de Diferenciación", xlab="Tiempo", ylab="Valores", col="2",lag=200)

#Cargamos la serie y las transformamos por logaritmos y luego se hace su diferencia

Glacial_varve=scan("varve.dat")
plot.ts(Glacial_varve, main="Sedimentos Glaciares", xlab="Tiempo", ylab="Valores", col="6")

acf(Glacial_varve,lag=200,main="Serie de No estacionaria-Datos Sedimentos Glaciares",col="2",xlab="Tiempo", ylab="Valores")
acf(diff(log(Glacial_varve),lag=1),lag=200,main="Transformación de lo datos de Sedimentos Glaciares",col="6",xlab="Tiempo", ylab="Valores")

Ejemplo_Rate-diferenciadasjpeg

Ejemplo_Rate-diferenciadas-Correlograma

Ejemplo_Sedimentos-GlciaresEjemplo_Sedimentos-Glciares-Correlograma

Ejemplo_Sedimentos-Glaciares_Transformados-Correlograma

Si bien no indiqué lo que es diferenciar una serie, basta notar que en el código cuando indico que se calcula la diferencial se reduce a considerar un retraso de la  seria contra el tiempo actual. Es decir, la diferencial en el tiempo t es y(t)-y(t-1).

Espero que las gráficas dejen claro cual es la idea de detectar un proceso estacionario o cuando es homogéneo de cierto grado.

Lo deseado en series de tiempo es construir un modelo tal que después de extraer la tendencia y la estacionalidad se puede tener certeza de que los residuos son muy parecidos a un ruido blanco gaussiano. Esto nos aseguro que obtuvimos toda la información posible desde los datos y que ya no se puede saber mucho más de la serie en los errores. También esto nos permite confiar en que el modelo es bueno y que sus predicción por ende deben de ser buenas.

Tomo un ejemplo de Avril Coghlan, el cual comenté en la entrada “series con tendencia y estacionalidad”. En esa entrada usé  el ejemplo para mostrar como usar los métodos Holt-Winters. En este caso lo uso para mostrar como se comportan los errores, los reviso gráficamente y aplicar la prueba Ljung-Box.

#Ejemplo análisis de residuos
#Extracción de datos y construcción del modelo Holt-Winters

souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987,1))
logsouvenirtimeseries <- log(souvenirtimeseries)
souvenirtimeseriesforecasts <- HoltWinters(logsouvenirtimeseries)
souvenirtimeseriesforecasts2 <- forecast.HoltWinters(souvenirtimeseriesforecasts, h=48)

#Revisión de la función de autocorrelación

acf(souvenirtimeseriesforecasts2$residuals, lag.max=20,main="Ejemplo de Correlograma de residuos de un modelo",xlab="Tiempo",ylab="Valores",col="6")

#Prueba Ljun-Box
Box.test(souvenirtimeseriesforecasts2$residuals, lag=20, type="Ljung-Box")

#La prueba muestra que se acepta la hipótesis nula, es decir, que son un ruido blanco
 Box-Ljung test

data: souvenirtimeseriesforecasts2$residuals
X-squared = 17.5304, df = 20, p-value = 0.6183

#Análisis gráfico

res=souvenirtimeseriesforecasts2$residuals
par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,30)
hist(res,prob=T,col="6")
lines(density(res),col="4")

Ejemplo_Correlograma_ResiduosEjemplo_análisis_gráfico_residuos

Observando el análisis gráfico y la prueba de hipótesis se concluye que el modelo no es malo y que los residuos se comportan como un ruido blanco gaussiano. Por lo cual se puede considerar que las predicción que se obtengan del modelo son buenas.

Modelos ARMA

La siguiente parte de la entrada trataré de mostrar con varios ejemplos aspectos de Medias Móviles (MA), Modelos Autorregresivos (AR), ARMA y SARMA.

Muchos aspectos teóricos solo los menciono y espero que los ejemplos ayuden a aprender el proceso de construcción de estos modelos. Evité poner ecuaciones en la medida de lo posible, pero algunas son obligadas debido a que es necesario para explicar algunos detalles.

Creo que comentaré algo más parecido a un recetario para identificar los procesos y solo menciono los nombres de los resultados teóricos que justifican los pasos, dejo algunas referencias o ligas para ver más detalles al respecto.

Suponemos que tenemos una serie de tiempo Y(t) y el operador de retrasos L() hace lo siguiente; tomando la serie en el tiempo t nos regresa el valor de la serie en el tiempo t-1. Es decir, L(Y(t))=Y(t-1).

Este operador por simple que parece nos permite construir teóricamente los modelos que comentaré y tener una notación simplificada. Se conoce como Lag Operator y pueden consultar sus propiedades básicas en las notas de Ruey Tsay.

Medias Moviles, MA(q). La idea de estos procesos o modelos es que dada nuestra serie Y(t) de algún experimento o fenómeno, sus valores son el resultado de sumas de un ruido blanco. Si nuestro ruido blanco se ve como e(t), entonces un modelo MA se ve como:

Y(t)=e(t)+B*e(t-1)+C*e(t-2)

donde se puede pensar e(t-1)=L(e(t)) y e(1-2)=L(L(e(t))), con el uso del operador de retrasos.

Sus propiedades son a grandes rasgos las siguientes:

  • Su media o esperanza es cero
  • La varianza del proceso en cualquier tiempo depende de los valores de los coeficientes B, C y de la varianza del ruido blanco.
  • Son procesos estacionarios.

La herramientas que siempre se usan para analizar series son la ACF y PACF, funcione de autocorrelación y función de autocorrelación parcial, respectivamente. La ACF de un MA(q) decrece rápidamente a cero y PACF decrece a cero si el procesos es invertible.

No menciono a detalle lo que significa que sea invertible, pero la idea es que si nuestro modelo de MA(q) tiene una estructura “manejable”, esto teóricamente nos permite deducir que los valores de la función  PACF se volverá cero cuando la revisemos para valores t muy grandes.

Una simulación de estos procesos se puede realizar con la librería forecast.

#Ejemplo Medias Móviles
#Cargamos la librería
libary(forecast)
n= 400
theta = c(1,-0.4,-0.4)

#Generamos la simulaicón

y.ma = arima.sim(list(order=c(0,0,2), ma=theta[2:3]), n=n, sd=sqrt(2.3))
layout(1:3)
ts.plot(y.ma,main="Modelo se MA(q)", xlab="Tiempo", ylab="Valore",col="6")
acf(y.ma,50, main="Medias Moviles")
pacf(y.ma,50, main="Medias Moviles")

Ejemplo_MA

Receta: revisar el comportamiento de la función de autocorrelación. La información de la existencia de un proceso de MA se concentra en ACF.

En el ejemplo tenemos un MA(2) y se muestra que las dos líneas que salen del la banda azul concuerdan con el número de parámetros del modelo y ACF converge a cero para t muy grande.

Modelos Autoregresivos, AR(q). Para estos modelos se considera que la información en Y(t) es debida a la información pasada de la serie y un error que se comporta como ruido blanco. Ejemplo:

Y(t)=A*Y(t-1)+B*Y(t-2)+…+Z*Y(t-p)+error(t)

De nueva el operador de retrasos L(), puede simplificar la notación.

Para asegurar que el proceso sea estacionario por lo menos en covarianza, se requiere cierta condición sobre los valores de A, B,..,Z, No doy detalles, pero pueden ser consultados en las notas de Ruey Tsay.

Otro aspecto teórico importante es el uso de la ecuación de Yule-Walker, la cual permite el cálculo de la función de Autocovarianza y su rol en la teoría y métodos de estimación es fundamental.

Receta: analizar la función ACFP.

La función ACF decae abruptamente indicando como las barras negras más largas la cantidad de parámetros del modelo y ACF decrece conforme t tiene a infinito. Esto caracteriza a los modelos autorregresivos.

 Unos ejemplos son los siguientes:

#Ejemplo 1 de modelo Autorregresivo
libary(forecast)
n= 400
theta = c(1,1.5,-0.9)
ts.sim <- arima.sim(list(order = c(2,0,0), ar = theta[2:3]), n = n)
layout(1:3)
ts.plot(ts.sim,main=expression("Modelo se AR(q) "*phi*"1=1.5 y "*phi*"2=-0.9"), xlab="Tiempo", ylab="Valore",col="6")
acf(ts.sim,50, main="Autoregresivo")
pacf(ts.sim,50, main="Autoregresivo")

#Ejemplo 2 de Autoregresivo
ts.sim <- arima.sim(list(order = c(1,0,0), ar = -.9), n = n)
layout(1:3)
ts.plot(ts.sim,main=expression("Modelo se AR(q) con "*phi*"=-0.9"), xlab="Tiempo", ylab="Valore",col="6")
acf(ts.sim,50, main="Autoregresivo")
pacf(ts.sim,50, main="Autoregresivo")

Ejemplo_AR1

Ejemplo_AR

Las gráficas muestran lo que se indicó, la función PACF indica el número de parámetros y ACF decae a cero.

Autorregresivos-medias móviles, ARMA(p,q). Estos modelos como el nombre indican consideran que la serie depende de la información pasada ( autorregresivo) y lo restante se modela por media de un proceso de medias móviles.

La idea es algo así:

Y(t)=AR(q)+Error(t)

donde Error(t)=MA(p)

Para asegurar que son estacionarios los procesos ARMA se requieren ciertas propiedades que son heredadas de su parte autorregresiva y de su parte de medias móviles. Al final de la entrada dejo algunas referencias.

Para aplicar la idea general es la siguiente, tenemos cierta serie Y(t) y modelamos su tendencia y su estacionalidad, la parte que resta del modelo son los residuos. Haciendo algunas pruebas estadísticas o analizando gráficamente los datos obtenemos que no son un ruido blanco o que podemos extraer más información. En esta última parte se procese a construir un modelo ARMA para los residuos.

El problema es elegir el modelo  ARMA(p,q) adecuado, no hay un método fijo, lo recomendable es hacer varios modelos candidatos. Para elegir uno de ellos es bueno considerar la parsimonia del modelo y apoyarse en los estadísticos AIC y BIC.

Esclareciendo lo anterior, trato de decir que no necesariamente el modelo con más parámetros es el mejor (parsimonia) y el comportamiento de AIC y BIC en la entrada “Series con tendencia y estacionalidad”  comenté como funcionan algunos estadísticos al momento de elegir entre varios modelos candidatos.

 Ejemplos de modelos ARMA

El primer ejemplo es solo para mostrar como se hace la construyen los modelos y al final se hace una simulación sobre los datos.

Dejo todo el código para que se replique el ejemplo en sus computadoras. Los datos son obtenidos desde la página.

#Ejemplo de modelos ARMA
#Análisis de Births
#Extracción de datos
library(forecast)
births <- scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
births.ts=ts(births[24:168],frequency=12,start=c(1948,1))

#Modelo 1-Holt-Winters
#Primero modelos de la serie

Mod1_Bir<-HoltWinters(births.ts)
Mod1_Bir_HW<-forecast.HoltWinters(Mod1_Bir)
res=ts(Mod1_Bir_HW$residual)

#Calculo de R cuadrada
y=births.ts
T = length(births.ts)
yest = fitted(Mod1_Bir_HW)
sse = sum((yest-y)^2)
ssr = sum((y-mean(y))^2)
R2 = 1-(sse/ssr)

#Análisis gráfico del modelo
plot.ts(births.ts, main="Modelo Holt-Winters",xlab="Tiempo",ylab="Valores")
lines(Mod1_Bir$fitted[,1],col=2)

#Análisis de los residuos
#Prueba Ljung-Box
Box.test(res, lag=20, type="Ljung-Box")
#Box-Ljung test

#data: res
#X-squared = 30.7366, df = 20, p-value = 0.05876
#La prueba indica que no se rechaza la hipótesis núla, que implica que los residuos sean un ruido blanco.
#Gráficas

par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,40)
hist(res,prob=T,col="6",nclass=20)
lines(density(res),col="4")

#Revisamos las autocorrelación y la autocorrelación parcial

par(mfrow=c(2,1))
acf(res,40)
pacf(res,40)

#Estimamos 3 modelos ARMA(2,2),ARMa(2,3),ARMA(3,2)
Mod1_res=arima(res,order=c(2,0,2))
Mod2_res=arima(res,order=c(2,0,3))
Mod3_res=arima(res,order=c(3,0,3))

#Por el criterio AIC elegimos ARMA(2,2)
Mod1_res$aic
Mod2_res$aic
Mod3_res$aic
#Revisamos sus comportamientos por medio de una prueba de diagnóstico

tsdiag(Mod1_res)
#########################################
#Modelo 2 quitamos tendencia y estacionalidad
births.ts2=ts(births[24:168])
t=1:length(births.ts)
sin1=sin(2*pi*t/12)
cos1=cos(2*pi*t/12)
#Modelo con tendencia lineas y estacionalidad
Mod1_TS=lm(births.ts2~t+sin1+cos1)

#Se revisa el modelo
#summary(Mod1_TS)

#Call:
#lm(formula = births.ts2 ~ t + sin1 + cos1)

#Residuals:
# Min 1Q Median 3Q Max 
#-3.4253 -0.6254 0.0933 0.6226 2.3010 

#Coefficients:
# Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 22.109601 0.173749 127.250 < 2e-16 ***
#t 0.044278 0.002066 21.434 < 2e-16 ***
#sin1 -1.050430 0.122438 -8.579 1.56e-14 ***
#cos1 -0.232017 0.121899 -1.903 0.059 . 
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 1.039 on 141 degrees of freedom
#Multiple R-squared: 0.7978, Adjusted R-squared: 0.7935 
#F-statistic: 185.5 on 3 and 141 DF, p-value: < 2.2e-16

#Gráfica del modelo
plot.ts(births.ts2)
lines(fitted(Mod1_TS),col=2)

#Revisión de los errores
res=Mod1_TS$residual
#Análisis de los residuos
#Se analiza gráficamentes y estadisticamente
#Prueba Ljung-Box
Box.test(res, lag=20, type="Ljung-Box")
#Box-Ljung test

#data: res
#X-squared = 30.7366, df = 20, p-value = 0.05876
#La prueba indica que no se rechaza la hipótesis núla, que implica que los residuos sean un ruido blanco.
#Gráficas

par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,40)
hist(res,prob=T,col="6",nclass=20)
lines(density(res),col="4")

#Revisamos las autocorrelación y la autocorrelación parcial

par(mfrow=c(2,1))
acf(res,40)
pacf(res,40)

#Estimamos 3 modelos ARMA(2,2),ARMa(2,3),ARMA(3,2)
Mod1_TS_res=arima(res,order=c(11,0,0))
Mod2_TS_res=arima(res,order=c(12,0,0))
Mod3_TS_res=arima(res,order=c(13,0,0))
Mod4_TS_res=arima(res,order=c(0,0,12))
Mod5_TS_res=arima(res,order=c(0,0,13))



#Por el criterio AIC elegimos ARMA(13,0)
Mod1_TS_res$aic
Mod2_TS_res$aic
Mod3_TS_res$aic
Mod4_TS_res$aic
Mod5_TS_res$aic

#Se elige el modelo Mod3_TS_res
#Revisamos sus comportamientos por medio de una prueba de diagnóstico

tsdiag(Mod3_TS_res)

#Gráfica de residuos
plot.ts(res)
lines(fitted(Mod3_TS_res))

#Cálculo de R cuadrada para los residuos

y=res
T = length(births.ts)
yest = fitted(Mod3_TS_res)
sse = sum((yest-y)^2)
ssr = sum((y-mean(y))^2)
R2 = 1-(sse/ssr)

#Simulación del modelo 
coeficientes=matrix(Mod3_TS_res$coef)
Mod_arma=arima.sim(list(order=c(13,0,0),ar=coeficientes[1:13,1]),n=145)

#Agregamos el valor del coeficientes constante
SimARMA=coeficientes[14,1]+Mod_arma
Simulacion=fitted(Mod1_TS)+SimARMA

#Graficas de los datos, la tendencia y estacionalidad y la simulación de modelos completo
plot.ts(births.ts2)
lines(fitted(Mod1_TS),col=2)
lines(Simulacion,col=6)

Algunas observaciones que se pueden hacer de los modelos anteriores, es que en el caso del modelo Holt-Winters no es necesario modelar los residuos por dos razones. La primera es que pasan las pruebas Ljung-Box, justo a penas pero la pasan y la segunda es que cuando se modelan los residuos y se elige uno de los modelos ARMA los pronósticos afectan poco  los pronósticos hechos por el método HW.

El segundo ajuste con tendencia y estacionalidad, si bien es malo permite jugar con los residuos. Al ser pocos datos no se puede esperar hacer un modelo más adecuado como ARMA(12,0,12) que es el que muestran los correlogramas. Entonces un ajuste medianamente bueno es ARMA(13,0,0) , pero se observa en la simulación que al comparar con los datos reales deja mucho que desear.

La intensión de poner todo el código anterior era mostrar lo mejor posible dos posibles procesos para elegir un modelo.

En el siguiente ejemplo se hacen tres tipos de modelos. Se analizan datos de ventas mensuales los cuales son extraídos desde la página de Robert H. Shumway.

#Ejemplo de construcción de modelos ARMA-HoltWinters y Análisis de frecuencias para modelo por medio de regresión lineal.
#Modelo de Sales
#Ejemplo 2 de modelo ARMA
#Extracción de lo datos
sales=scan("sales.dat")
plot.ts(sales,main="Ventas", ylab="Valor", xlab="Tiempo")
plot.ts(diff(sales,1),ylab="Valores diferenciados",xlab="Tiempo", main="Datos de Ventas diferenciados")
dif.sales=diff(sales,1)
#Exploración de las frecuencias relevantes
#Jugando con el espectro de los datos
P=abs(2*fft(dif.sales)/n)^2
f=0:75/150
plot(f,P[1:76],type="o",main="Espectro de Frecuencias",xlab="Frecuencias",ylab="Valores")
abline(v=0.012,col=2)
abline(v=0.04,col=2)
abline(v=0.06,col=2)
abline(v=0.093,col=2)
abline(v=0.26,col=2)

#Onteniendo las frecuencias se pueden obtener los periodos
t=1:length(dif.sales)
sin1=sin(2*pi*t/77)
cos1=cos(2*pi*t/77)
sin2=sin(2*pi*t/25)
cos2=cos(2*pi*t/25)
sin3=sin(2*pi*t/16.5)
cos3=cos(2*pi*t/16.5)
sin4=sin(2*pi*t/10.75)
cos4=cos(2*pi*t/10.75)
sin5=sin(2*pi*t/3.84)
cos5=cos(2*pi*t/3.84)
#Modelos con distintos periodos
#dif.sales=diff(sales,1)
Mod1_TS=lm(dif.sales~t+sin1+cos1)
Mod2_TS=lm(dif.sales~t+sin1+cos1+sin2+cos2)
Mod3_TS=lm(dif.sales~t+sin1+cos1+sin2+cos2+sin3+cos3)
Mod4_TS=lm(dif.sales~t+sin1+cos1+sin2+cos2+sin3+cos3+sin4+cos4)
Mod5_TS=lm(dif.sales~t+sin1+cos1+sin2+cos2+sin3+cos3+sin4+cos4+sin5+cos5)
#Graficas de ajustes
plot.ts(diff(sales,1),main="Ejemplo de Ajustes con tendencia")
lines(fitted(Mod1_TS),col=2)
lines(fitted(Mod2_TS),col=6)
lines(fitted(Mod3_TS),col=8)
lines(fitted(Mod4_TS),col=4)
lines(fitted(Mod5_TS),col=3)
#Se selecciona el modelo Mod5_TS
res=Mod5_TS$residual
#Análisis de los residuos
#Se analiza gráficamentes y estadisticamente
#Prueba Ljung-Box
Box.test(res, lag=20, type="Ljung-Box")
#La prueba indica que no se rechaza la hipótesis núla, que implica que los residuos sean un ruido blanco.
#Gráficas

par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,40)
hist(res,prob=T,col="6",nclass=20)
lines(density(res),col="4")

#Revisamos las autocorrelación y la autocorrelación parcial

par(mfrow=c(2,1))
acf(res,40,col=2)
pacf(res,40,col=2)

#Jugando con el método Holt-Winters
Mod1_HW<-HoltWinters(ts(dif.sales,frequency=3))
Mod1_sales_HW<-forecast.HoltWinters(Mod1_HW)

plot.ts(ts(dif.sales,frequency=3), main="Modelo Holt-Winters",xlab="Tiempo",ylab="Valores")
lines(fitted(Mod1_sales_HW),col=2)

res=Mod1_sales_HW$residual
#Análisis de los residuos
#Se analiza gráficamentes y estadisticamente
#Prueba Ljung-Box
Box.test(res, lag=20, type="Ljung-Box")
#La prueba indica que no se rechaza la hipótesis núla, que implica que los residuos sean un ruido blanco.
#Gráficas

par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,40)
hist(res,prob=T,col="6",nclass=20)
lines(density(res),col="4")

# Modelos ARMA sobre las ventas diferenciadas
#Exploramos los posibles modelos
par(mfrow=c(2,1))
acf(dif.sales,40,col=2)
pacf(dif.sales,40,col=2)

#Calculamos varios modelos ARMA
Mod1_ARMA=arima(dif.sales,order=c(1,0,1))
Mod2_ARMA=arima(dif.sales,order=c(1,0,2))
Mod3_ARMA=arima(dif.sales,order=c(1,0,3))
Mod4_ARMA=arima(dif.sales,order=c(1,0,4))
Mod5_ARMA=arima(dif.sales,order=c(2,0,1))
Mod6_ARMA=arima(dif.sales,order=c(2,0,2))
Mod7_ARMA=arima(dif.sales,order=c(2,0,3))
Mod8_ARMA=arima(dif.sales,order=c(2,0,4))

#Elección del modelo
#Se elige el modelo Mod1_ARMA
Mod1_ARMA$aic
Mod2_ARMA$aic
Mod3_ARMA$aic
Mod4_ARMA$aic
Mod5_ARMA$aic
Mod6_ARMA$aic
Mod7_ARMA$aic
Mod8_ARMA$aic

plot.ts(diff(sales,1),ylab="Valores diferenciados",xlab="Tiempo", main="Datos de Ventas diferenciados")
lines(fitted(Mod1_ARMA),col=2)

res=Mod1_ARMA$residual
#Análisis de los residuos
#Se analiza gráficamentes y estadisticamente
#Prueba Ljung-Box
Box.test(res, lag=20, type="Ljung-Box")

#La prueba indica que no se rechaza la hipótesis núla, que implica que los residuos sean un ruido blanco.
#Gráficas

par(mfrow=c(2,2))
plot.ts(res)
abline(h=0,lty=2,col="2")
qqnorm(res)
qqline(res,col="2")
acf(res,40)
hist(res,prob=T,col="6",nclass=20)
lines(density(res),col="4")

Espectro_de_Frecuencias_Ventas_dif

Ejemplos_Ajustes_por_Periodos

 En estos últimos modelos se hace uso de cosas que comente en otras entradas, ejemplo; analizar el espectro de frecuencias para construir un modelo y diferenciar la serie original para analizarla. Algo que no menciono a detalle pero que aparece en los modelos es que al diferenciar y ajustar por un modelo ARMA, lo que realmente se realiza es un modelo ARIMA. Sobre esos modelos comento en otra ocasión.

 Referencias:

1.-http://anson.ucdavis.edu/~shumway/

2.-https://www.otexts.org/fpp/

3.-http://www.statsoft.com/Textbook/Time-Series-Analysis

4.-http://www.medellin.unal.edu.co/~ndgirald/members.htm

5.-http://pandas.pydata.org/pandas-docs/stable/timeseries.html

6.-Analysis of financial time series

7.-Nonlinear time series analysis

8.-Econometric Models Economic Forecasts

Anuncios

Series con tendencia y estacionalidad.

En esta entrada trataré de ejemplificar los métodos que se usan cuando se tienen una serie con tendencia y/o  estacionalidad.

Series con tendencia.

Hablar de tendencia es más o menos evidente cuando uno ve la gráfica de la serie de tiempo. Se observa, o una posible recta o una curva y lo que uno hace es ajustar los datos o transformar los valores de la serie, ejemplo por logaritmos; para ver algún posible ajuste.

Pero detrás de la idea de un ajuste, está la idea de hacer una regresión a los datos, ya sea lineal o no lineal, es el caballito de batalla al momento de hacer algún modelo o alguna predicción. Entonces las herramientas son estándar y bien conocidas.  Para ejemplificar tomo dos ejemplos de Robert H. Shumway.

#Ejemplo 1
#Se extraen los datos desde la página R. Shumway

temperatura_global<-scan("http://anson.ucdavis.edu/~shumway/globtemp.dat")
x=temperatura_global[45:142]
t=1900:1997
#Estimamos la regresión o tendencia
fit=lm(x~t)
#Revisamos el cálculo de la tendencia
summary(fit)

Call:
lm(formula = x ~ t)

Residuals:
 Min        1Q     Median   3Q     Max 
-0.30352 -0.09671 0.01132 0.08289 0.33519 

Coefficients:
           Estimate   Std. Error t value Pr(>|t|) 
(Intercept) -1.219e+01 9.032e-01 -13.49 <2e-16 ***
t             6.209e-03 4.635e-04 13.40 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1298 on 96 degrees of freedom
Multiple R-squared: 0.6515, Adjusted R-squared: 0.6479 
F-statistic: 179.5 on 1 and 96 DF, p-value: < 2.2e-16

plot(t,x,type="o",xlab="años",ylab="Temperatura",col="2",main="Ejemplo de cálculo de Tendencia")

abline(fit,col="4")

 Ejeplo1_de_Tendencia

En este ejemplo se ve como al hacer el ajuste la recta parece ser un buen modelo de tendencia, y el resumen de estadístico del ajuste indica que se explica el 64% de variabilidad de la serie y se puede considerar como valido el ajuste al ver el p-Valor. Sobre estos estadísticos comento más adelante.

En el segundo ejemplo considero 4 posibles modelos entre tres series de tiempo. La serie que tratamos de explicar o modelar indica la media de muertes diarias por enfermedades cardiovasculares y trata de explicar su comportamiento mediante los datos de la contaminación y  temperatura. Los datos son correspondientes a un periodo de años desde 1990 a 1999.

Este modelo bien puede considerarse como parte de una herramienta de exploración , pero permite ver dos cosas más. La primera, mostrar como sugerir varios modelos y la segunda evaluar los modelos revisando sus estadísticos. Esto con la finalidad de elegir uno de ellos como el “mejor”.

#Ejemplo 2
#Se Cargan los datos
mortalidad<-scan("http://anson.ucdavis.edu/~shumway/cmort.dat")
Temp<-scan("http://anson.ucdavis.edu/~shumway/temp.dat")
Pollution<-scan("http://anson.ucdavis.edu/~shumway/part.dat")

t=1:lenght(mortalidad)

#Modelo 1
#Se ajusta la serie contra el tiempo
#Mortalidad=A+B*t+error
fit1<-lm(mortalidad~t)
summary(fit1)

Call:
lm(formula = mortalidad ~ t)

Residuals:
   Min    1Q    Median 3Q    Max 
-16.445 -6.670 -1.366 5.505 40.107 

Coefficients:
          Estimate   Std. Error  t value Pr(>|t|) 
(Intercept) 96.651347 0.790318   122.29 <2e-16 ***
       t    -0.031247 0.002691   -11.61 <2e-16 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.893 on 506 degrees of freedom
Multiple R-squared: 0.2104, Adjusted R-squared: 0.2089 
F-statistic: 134.9 on 1 and 506 DF, p-value: < 2.2e-16

#Modelo 2
#Mortalidad=A+B*t+c*(T-mean(T))+error
#T=Temperatura
mean(Temp)
74.26041
T=mean(Temp)
Temp2=Temp-T
fit2=lm(mortalidad~t+Temp2)
summary(fit2)

Call:
lm(formula = mortalidad ~ t + Temp2)

Residuals:
   Min    1Q    Median  3Q   Max 
-16.846 -5.330 -1.207 4.701 33.306 

Coefficients:
            Estimate  Std. Error t value Pr(>|t|) 
(Intercept) 96.22547 0.70183   137.11 <2e-16 ***
       t    -0.02957 0.00239   -12.37 <2e-16 ***
      Temp2 -0.45792 0.03893   -11.76 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.887 on 505 degrees of freedom
Multiple R-squared: 0.3802, Adjusted R-squared: 0.3778 
F-statistic: 154.9 on 2 and 505 DF, p-value: < 2.2e-16

#Modelo 3
#Mortalidad=A+B*t+c*(T-mean(T))+d*(T-mean(t))^2+error
T2=Temp^2
fit3=lm(mortalidad~t+Temp2+T2)
summary(fit3)
lm(formula = mortalidad ~ t + Temp2 + T2)

Residuals:
  Min     1Q    Median  3Q   Max 
-17.464 -4.858 -0.945 4.511 34.939 

Coefficients:
             Estimate  Std. Error t value Pr(>|t|) 
(Intercept) 93.918196 0.725174   129.511 < 2e-16 ***
      t     -0.028738 0.002261   -12.710 < 2e-16 ***
      Temp2 -0.480793 0.036895   -13.031 < 2e-16 ***
         T2  0.025830 0.003287    7.858 2.38e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.452 on 504 degrees of freedom
Multiple R-squared: 0.4479, Adjusted R-squared: 0.4446 
F-statistic: 136.3 on 3 and 504 DF, p-value: < 2.2e-16

#Modelo 4
#Mortalidad=A+B*t+C*(T-mean(T))+D*(T-mean(t))^2+Poll+Error

fit4=lm(mortalidad~t+Temp2+Temp2^2+Pollution)
summary(fit4)
Call:
lm(formula = mortalidad ~ t + Temp2 + T2 + Pollution)

Residuals:
   Min     1Q     Median  3Q     Max 
-19.0760 -4.2153 -0.4878 3.7435 29.2448 

Coefficients:
             Estimate Std. Error t value Pr(>|t|) 
(Intercept)   81.592238  1.102148 74.03  < 2e-16 ***
         t    -0.026844  0.001942 -13.82 < 2e-16 ***
        Temp2 -0.472469 0.031622  -14.94 < 2e-16 ***
         T2    0.022588 0.002827   7.99 9.26e-15 ***
     Pollution 0.255350 0.018857  13.54 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.385 on 503 degrees of freedom
Multiple R-squared: 0.5954, Adjusted R-squared: 0.5922 
F-statistic: 185 on 4 and 503 DF, p-value: < 2.2e-16

De estos 4 modelos se tienen suficientes datos para revisar los estadísticos  importantes sobre una regresión.

1.-R^2  o R cuadrada,  es un estadístico descriptivo y que su valor representa el porcentaje de variabilidad descrita por las variables explicativas. En general se toma la R cuadra ajustada para analizar el resultado de la regresión. En los modelos se observa que el Modelo 1 describe solo el 20%, el Modelo 2 el 37%, el Modelo 3 el 44% y el modelo 4 el 59%.

La receta para aceptar o rechazar un modelo revisando R cuadrada es más o menos así, si su valor es menor a 40% es un modelo de regular a malo, si es arriba de 60% es bueno y entre 40% y 60% es recomendable revisar el modelo para ver si agregando una nueva variable o transformando alguna se mejora el porcentaje.

Al ser un estadístico descriptivo un bajo valor implica no directamente rechazar el modelo, sino explorar los datos y validar si se puede mejorar agregando nuevas variables o transformando los datos.

2.-El p-Valores, indica que en caso de ser menor de 0.05 el suponer que los valores de los coeficientes de las variables sean cero se debe de rechazar. Es decir; la hipótesis nula de que el valor del coeficiente es cero se rechaza en todos los modelos anteriores.

3.-El F-estadístico, indica que se realiza una prueba de hipótesis suponiendo que todos los coeficientes son nulos. En R project se hace una el cálculo del p-Valor  de este estadístico, observamos que en todos los casos se rechaza la hipótesis nula, así no todos los coeficientes son distintos de cero.

4.-Un estadístico que no se calcula directamente de la regresión en R project es AIC (Akaike’s Information Criterion). La idea de este es elegir el modelo que minimiza su valor. Se calculan sus valores y se elige un modelo. Otro estadística que funciona similar AIC es BIC( Bayesian Information Criterion) o SIC (Schwarz’s Information Criterion), no lo calculo en este ejemplo pero la función en R es BIC.

#Calculo de AIC
 AIC(fit1)
3665.899
AIC(fit2)
3544.889
AIC(fit3)
3488.176
AIC(fit4)
3332.282

Se observa que el modelo 4 es el menor valor de AIC  al compararlo con respecto a los otros. Entonces, se elige el modelo 4 como el mejor modelo al revisar los estadísticos.

 Ejemplo2_regresion_lineal

 Los dos ejemplos anteriores muestran como estimar una tendencia y cuales son los estadísticos usuales para definir entre posibles modelos de tendencia el más adecuado.

El tipo de tendencias que se suelen modelar son: lineal, cuadrática, cúbica, exponencial, Log-Lineal, polinómica.

El polinómico es una generalización de cuadrático y cúbico, lo que puede suceder es que sobre estime la tendencia ya que nada restringe elegir un polinomio de grado 4 o 20. En varios casos cuando se usan este tipo de modelos se usan la regularización y la validación cruzada para determinar el mejor modelos entre los polinómicos.

Los modelos consideran la estimación de la tendencia sobre todos los datos, es como extraer la tendencia global de la serie. Esto es restrictivo o muy rígido al modelar la series, en algunos casos prácticos es suficiente este tipo de estimación de la tendencia. Para otras situaciones se cuenta con métodos no-paramétricos que hacen el cálculo de la tendencia de manera loca. Ejemplo, puede suceder que localmente los datos muestren tendencias lineales pero donde los parámetros de la recta a(t)+b(t)*T; van variando con el tiempo.

Tomando los datos de Ruey S. Tsay del indicador Standar&Poor’s de Enero de 1950 hasta Abril 2008  estimo 4 tipo de tendencias y comparo los estadísticos para elegir alguna como el mejor modelo de la tendencia.

#Ejemplo 3
#Extraemos los datos
sp<-read.table("d-sp55008.txt",head=T)
head(sp)
sp=sp[,7]
#Definimos la variable tiempo
t=1:length(sp)
t2=t^2
#Hacemos las estimaciones de los modelos
#Cuadrático
fit.cuadratico=lm(sp~t+t2)
#Calculamos el logaritmo de la serie
log.sp=log(sp)
#Estimamos este ajuste como auxiliar
fit.Log.Lin=lm(log.sp~t)

#Usamos ciertos datos para estimar el mod. exponecial
b_0=fit.Log.Lin$coefficient[1]
b_1=fit.Log.Lin$coefficient[2]
Datos=data.frame(sp,t)
#Calculamos por regresión no lineal el modelo exponencial
mod.fit.exp=nls(sp~exp(beta0+beta1*t),data=Datos,start=list(beta0=beta_0,beta1=beta_1))
#Log.Lineal
Log.Lin.Mod=exp(a_0+a_1*t)
#Lineal 
fit.lineal=lm(sp~t)

#Hacemos la gráfica para ver los ajustes
plot(ts(sp),xlab="Desde 1950 hasta 2008",ylab="Valores de S&P",main="Ejemplo de calculo de Tendencias")
lines(predict(fit.lineal),col="2")
lines(predict(mod.fit.exp),col="3")
lines(Log.Lin.Mod,col="4")
lines(fit.cuadratico$fitted,col="6")
legend("topleft",c('Lineal','Exponencial','Log-Lineal','Cuadrática'),col=c(2,3,4,6),pch=1,text.col=c(2,3,4,6))


 Ejemplo_Ajuste_Tendencias

Lo que ahora hago es calcular las medidas estadísticas para cada modelo, usamos una función que concentra todos los estadísticos.

#Elección del modelo
#Función con todas las medidas
medidas = function(m,y,k){
# y = serie, m = modelo, k = numero parametros
T = length(y)
yest = fitted(m)
sse = sum((yest-y)^2)
ssr = sum((y-mean(y))^2)
mse = sse/(T-k)
R2 = 1-(sse/ssr)
Ra2 = 1-((T-1)*(1-R2)/(T-k))
aic = log((T-k)*exp(2*k/T)*mse/T)
bic = log(T^(k/T)*(T-k)*mse/T)
M = c(Ra2, mse, aic, bic)
names(M) = c("R2-ajus","MSE","logAIC","logBIC")
return(M)}

#Elección del modelo
m.lin=medidas(fit.lineal,sp,2)
m.exp=medidas(mod.fit.exp,sp,2)
m.cuad=medidas(fit.cuadratico,sp,3)
m.LogLin=medidas(fit.Log.Lin,log(sp),2)
M=cbind(m.lin,m.exp,m.cuad,m.LogLin)
> M
           m.lin      m.exp        m.cuad    m.LogLin
R2-ajus 6.911007e-01 9.126325e-01 9.087651e-01 0.94441124
    MSE 5.889685e+04 1.665809e+04 1.739547e+04 0.08440446
 logAIC 1.098368e+01 9.720788e+00 9.764170e+00 -2.47199860
 logBIC 1.098472e+01 9.721823e+00 9.765723e+00 -2.47096286

 Se hace notar inmediatamente que el modelo que resulta tener mejores estadísticos es el Log-Lineal. Lo cual desde la gráfica se visualizaba un poco.

Una ves elegido el modelo para la tendencia es necesario revisar los residuos o lo que se identifica como su parte irregular. Solo pongo las gráficas estándar para explorar su comportamiento.

#Revisión de residuos
par(mfrow=c(2,2))
plot.ts(res.LogLin)
abline(h=0,lty=2,col="2")
qqnorm(res.LogLin)
qqline(res.LogLin,col="2")
acf(res.LogLin,30)
hist(res.LogLin,prob=T)
lines(density(res.LogLin),col="4")

 Gráficas_de_los_errores

En las gráficas de los residuos se observa que el ajuste no es muy adecuado. En la última parte de esta entrada hablo al respecto de condiciones necesarias de los residuos para considerar un buen ajuste.

Alisado de una serie 

Un par de métodos para alisar series y para hacer pronósticos desde un modelo determinista son: medias móviles y medias móviles ponderadas exponencialmente.

En caso de usar esos métodos para hacer pronósticos la idea es considerar que la media de una cantidad de datos de la serie es un buen pronóstico , esto en el caso de las medias móviles. Para el caso de las medias móviles ponderadas exponencialmente se toma la misma idea, que un promedio de datos de la serie pueden dar un buen pronóstico, pero comparado con las medias móviles se considera que los últimos datos de la serie tienen mayo peso al momento de ofrecer un pronóstico.

La idea de alisar una serie es visualizar con claridad la tendencia y patrones críticos de la serie. Se trata hasta cierto punto de resaltar el comportamiento más relevante de la serie sin formular de forma explicita un modelo  o hipótesis.

Para hacer unos ejemplos de suavizado tomo datos desde la página de Ruey S.Tsay. Estos datos corresponden a las tasas de interés mensual de EEUU.

#Alisado de una serie

#Medias Móviles
#Extraemos los datos
datos=read.table("m-gs10.txt",head=T)
datos1=datos$value[1:300]
#Alisamos MA centrados
alisado1=filter(datos1,rep(1,12)/12,side=2)
alisado2=filter(datos1,rep(1,24)/24,side=2)

#Hacemos su gráfica para comparar
plot.ts(datos1,main="Serie de tasas de interes",xlab="Tiempo",ylab="Valores")
lines(alisado1,col="2")
lines(alisado2,col="4")

#Alisado MA no centrado
alisado3=filter(datos1,rep(1,12)/12,side=1)

#Alisado exponencial 
#Usamos HoltWinters omitiendo el valor de beta y gamma

alisado5=HoltWinters(datos1,beta=FALSE,gamma=FALSE)
alisado6=HoltWinters(datos1,alpha=0.05,beta=FALSE,gamma=FALSE)
plot.ts(datos1,main="Ejemplo de alisados para las tasas de interes",xlab="Tiempo", ylab="Valores")
lines(alisado5$fitted[,1],col="2")
lines(alisado6$fitted[,1],col="4")
legend("topleft",c('Datos','Alisado sin alpha dado','Alisado con alpha 0.05'),col=c(1,2,4),pch=1,text.col=c(1,2,4))

#Alisado exponencial incorporando cambios promedios en la tendencia

alisado8=HoltWinters(datos1,gamma=FALSE)
alisado9=HoltWinters(datos1,alpha=0.2,beta=.45,gamma=FALSE)
plot.ts(datos1,main="Ejemplo de alisados para las tasas de interes",xlab="Tiempo", ylab="Valores")
lines(alisado8$fitted[,1],col="4")
lines(alisado9$fitted[,1],col="6")
legend("topleft",c('Datos','Alisado sin alpha y beta datos','Alisado con alpha 0.2 y beta=0.45'),col=c(1,4,6),pch=1,text.col=c(1,4,6))

 En la primera gráfica se observa la estimación de dos medias móviles centradas, para 12 y 24 datos.

Alisado_ma_1

En la siguiente gráfica se observa como se alisa la serie con medias móviles centradas y no centradas para la misma cantidad de datos, 12 observaciones de la serie.

Alisado_ma_2

La tercera gráfica muestra el método de alisado exponencial, existen varias librerías para estimar esto, pero se puede usar la función HoltWinters que es parte de las funciones base de R project. Cuando no le cargamos el valor de alpha hace la estimación de su valor estimando las soluciones recurrentes del modelo de suavizado exponencial.

Alisado_mape1

La cuarta gráfica muestra el alisado con el método Holt, lo cual también es un método de alisado exponencial pero considera cambios promedios en la tendencia de la serie. Para su estimación en R se usa la función HoltWinters omitiendo solo el valor de gamma. Cuando no definen los valores de alpha y beta, la función hace el calculo de los valores haciendo uso de las ecuaciones recurrentes del método.

La regla para determinar el uso de alpha en los dos últimos métodos que son de suavizado exponencial es, una alpha cercana a cero hace un alisado con menos variaciones o curvas. Con alpha cercana a 1 indica que el peso de los últimos valores de la serie tienen mayor peso en el alisado.

 Estos métodos de alisado se pueden usar para hacer predicciones a corto plazo, como ejemplo calculo los siguientes 12 valores de la serie.

#Ejemplo de pronóstico de 12 valores futuros
#Para hacer los pronósticos de MA uso la siguiente ffunción

Predict_MA<-function(serie,num){ 
 #Funciona bien para más de 10 valores a pronósticar 
 L=length(datos1)
 x=datos1[(L-num):L]
 v=x
 predic=matrix(nrow=1,ncol=1) 
 for( i in (1:num)){
 y=(sum(v)/num)
 v=c(x[(i+1):num],v[(num-i):11],y)
 predic=rbind(predic,y)
 }
 return(predic[2:(num+1)])}


#Gráficas de pronósticos

plot.ts(datos1,xlim=c(1,315),xlab="tiempo",main="Ejemplo de Pronóstico", ylab="Valores")
lines(predict(alisado5,n.ahead=12),col="6")
lines(alisado5$fitted[,1],col="6")
lines(c(alisado3,V),col="4")
legend("topleft",c('Datos','Alisado Exponencial','Alisado Media Móvil'),col=c(1,6,4),pch=1,text.col=c(1,6,4))

 Se tiene la siguiente gráfica como resultado del código anterior.

Ejemplo_Pronostico_MA_Exp_12v

La librería que se usa con mayor frecuencia para hacer pronósticos por medio de las series de tiempo es forecast. En otra entrada explico como usarla y sobre todo, como emplearla para modelos estocásticos de series de tiempo.

También para realizar el suavizado exponencial existe la librería fTrading y otra opción es la librería signal.

 Estacionalidad de la serie de tiempo

Cuando se piensa en la serie de tiempo formada por sus componente no observables, la estacionalidad indica el movimiento cíclico en cierto periodo de tiempo el cual se repite sistemáticamente. Ejemplo, si la estacionalidad es una función del tiempo S(t) y se repite cada cierto tiempo m, entonces S(m+t)=S(t), por lo cual es una función periódica.

El periodo en el cual se repite el patrón estacional depende del origen de los datos, puede ser anual, trimestral, mensual, semanal o diario.

Para modelar o estimar la parte estacional de una serie se puede hacer un modelo considerando variables indicadoras o funciones trigonométricas (seno y coseno). En el caso de las variables indicadoras, si la estacionalidad se muestra cada semana y la serie representa datos de todo un año tendríamos 52 variables indicadoras, lo cual es conveniente modelar la serie con funciones trigonométricas en lugar de indicadoras.

 Un ejemplo de una serie con estacionalidad y que puede ser modelada con componentes aditivos, es debida a la demanda diaria de energía. Hacemos dos modelos, uno con variables indicadoras o dummy y el segundo con funciones trigonométricas. Uso la librería forecast para estimar variables generar las variables indicadoras y las trigonométricas. Tomo unos datos de Norma Giraldo Gómez.

#Ejemplos de Cálculo de la estacionalidad
#Extraemos los datos y solo usamos la mitad de ellos
demanda<-read.table("http://www.medellin.unal.edu.co/~ndgirald/Datos/Datos%20curso%20Series%20II/demanda_energia_18.dat",header=TRUE)
dem=demanda[1:1461]
t=1:length(dem)
dem=ts(dem,frequency=7)

#Cargamos la librería forecast
library(forecast)

#Modelo 1, combinamos indicadoras con trigonométricas para detectar la estacionalidad anula y semanal

St=seasonaldummy(dem)
X1=cos(2*pi*t/365)
X2=sin(2*pi*t/365)
mod1=lm(dem~t+St+X1+X2)

#Modelo 2, generamos una matriz de funciones trigonométricas.

Trig=fourier(dem,2)
mod2=lm(dem~t+Trig)

#Graficamos los resultados

plot.ts(dem,main="Ejemplo de serie con estacionalidad",xlab="Tiempo", ylab="Valore de la demanda de energía")
lines(mod1$fitted,col="2")
lines(mod2$fitted,col="4")

Ejemplo_Estacionalidad_1

 En los dos modelos se consideró que la tendencia era lineal y se modelo de dos formas la estacionalidad. Pese a que la gráfica no muestre con mucho detalle el comportamiento de los modelos se puede revisar el resumen de los datos y se observará que el primero explica el 63% de variabilidad y el segundo 57%.

Debido a que se modela la tendencia y la estacionalidad, estos modelos se llaman estructurales. Lo que se hace para validar si los modelos son adecuados es analizar el comportamiento de los coeficientes, debido a que puede suceder que los parámetros estimados cambien conforme al tiempo.

Si los parámetros muestran ese comportamiento no son un buen modelo para fiarse al hacer predicciones. Dos herramientas para revisar el comportamiento de los coeficientes es por medio de sus estimación de manera recursiva o mediante la prueba CUSUM.

#Ejemplo de análisis de los coeficientes
#Primero se hacen estimaciones del modelo con funciones trigonométricas recursivamente

s=4
k = 2+3+10
T=length(dem)
n = T-k
parm = mat.or.vec(n,(2+s))
for(j in (1:n)){
y = dem[1:(k+j)]
t1 = t[1:(k+j)]
It1 = Trig[1:(k+j),]
mod = lm(y~ t1+It1)
parm[j,] = t(mod$coefficient)
}
colnames(parm)=c("beta.0","beta.1",
"delta.1","delta.2","delta.3","delta.4")
#Se grafica el comportamiento de cada uno de los coeficientes
plot.ts(parm)

#Ejemplo de prueba CUSUM
#Se usa la librería strucchange

library(strucchange)
#Se calculan los residuos
rrc = recresid(dem ˜ t + Trig)
#Se hacen pruebas gráficas por dos métodos de estimación de los residuos

prueba.cusum = efp( dem~t+Trig, type = "RecCUSUM")
plot(prueba.cusum1)
prueba.cusum2 = efp(dem~t+Trig, type = "OLSCUSUM")
plot(prueba.cusum2)

 Las dos gráficas que se obtienen son las siguientes:

Estabilidad1_recursivo

Ejem_CUSUM

Lo que se observa es que el modelo por funciones trigonométricas no tienen coeficientes que sean constante conforme varia el tiempo. Es decir, la estabilidad del modelo es rechazada.

 Otro tipo de técnicas que se usa para analizar la series cuando tienen estacionalidad son los métodos de ajuste estacional, que se usan para calcular el índice de estacionalidad y una vez calculados desestacionalizar la serie. Sobre esto no comento nada, se pueden buscar las técnicas y sobre todo son usadas en econometría.

El último método es llamada Holt-Winters, este método estima recursivamente los parámetros del modelo considerando que tienen una tendencia lineal. Tanto se puede usar para alisar la serie con estacionalidad como para realizar pronósticos.

En R project estimar el modelos Holt-Winters es sencillo ya que tienen una función por default, pero también se pueden hacer estimaciones o pronósticos con la librería forecast.

Muestro dos ejemplos, uno de la serie de demanda de energía y otro de ventas de una tienda. El segundo ejemplo se puede consultar en las notas de Avril Coghlan. Pongo este último ejemplo por que creo que visualmente es más claro como se estima y la facilidad con la cual se puede hacer un buen modelo con la técnica Holt-Winters.

#Ejemplo Holt-Winters
dem.hw<-HoltWinters(dem)
plot.forecast(forecast.HoltWinters(dem.hw,h=48))

#Ejemplo 2
#Se debe de tener instalada la libreria forecast

souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987,1))
logsouvenirtimeseries <- log(souvenirtimeseries)
plot.ts(logsouvenirtimeseries)
souvenirtimeseriesforecasts <- HoltWinters(logsouvenirtimeseries)
plot(souvenirtimeseriesforecasts)
souvenirtimeseriesforecasts2 <- forecast.HoltWinters(souvenirtimeseriesforecasts, h=48)
plot.forecast(souvenirtimeseriesforecasts2)

 Ejemplo_HoltWinters2

Ejemplo_HoltWinters3

 Espero que la entrada de una idea general de las técnicas, algunas cosas que no menciono son las técnicas Loess, técnicas de Kernel y vecinos cercanos. Son técnicas útiles para exploración de la tendencia de la serie o el comportamiento en general.

Todo el material puede ser consultado en alguna buena referencia  para conocer la teoría al respecto a estos temas, me reservo comentar alguna referencia ya que pese a que son las mismas técnicas ver un libro enfocado a finanzas es distinto a ver la perspectiva de uno de geofísica.

Referencias:

1.-http://anson.ucdavis.edu/~shumway/

2.-https://www.otexts.org/fpp/

3.-http://www.statsoft.com/Textbook/Time-Series-Analysis

4.-http://www.medellin.unal.edu.co/~ndgirald/members.htm

5.-http://pandas.pydata.org/pandas-docs/stable/timeseries.html

6.-Analysis of financial time series

7.-Nonlinear time series analysis

8.-Econometric Models Economic Forecasts

Exploración de una Serie de Tiempo

Sobre Exploración de una serie de tiempo

Para hacer la exploración de la serie se hace uso de las herramientas básica de exploración de datos: medidas de localización, medidas de dispersión y de forma.

Para las medidas de localización las funciones que se requieren son las usuales mean,quantile,summary. Para medidas de dispersión se usan  range,IQR,var,sd y una opción para las medidas de forma es usar la librería e1071 donde podemos calcular las funciones kurtosis,skewness.

En resumen estas funciones nos dan información de la densidad, de la distribución y de como se comportan los datos en general. Siempre la mejor herramienta para tener una idea rápida sobre los datos es alguna gráficas.

#Analizamos los datos
glotemp<-scan("glotemp.dat")
hist(glotemp,col="2",main="Histograma",xlab="Temperatura global",ylab="Densidad",prob=T)
lines(density(globtemp),lwd=2)

qqnorm(glotemp,main="QQ-Norm",xlab="Cuantiles teóricos",ylab="Cuantiles muestrales",col="2")

 Histograma_densidad_temp-global

QQ-norm_datos_temperaturas_globales

Esto no es diferente a la exploración de datos en general, pero como al analizar una serie buscamos si existe algún comportamiento solo desde la serie misma, resulta importante explorar como se comporta la serie al graficar contra si con tiempo diferentes. Es decir, analizamos la serie completa contra la serie retardada un tiempo para ver si existe algún patrón de comportamiento, pero podemos hacer este análisis contra tiempos retrasados mayores.

El siguiente ejemplo muestra lo que trato de mostrar en el párrafo anterior.

#Scatterplot de serie contra retrasos
soi<-scan("soi.dat")
lag.plot(soi,lags=12,layout=c(3,4),diag=F,colo="2")

 Los datos soi corresponden al Southern Oscillation Index y los datos se pueden descargar desde la página de Robert Shumway.

También puede suceder que sospechamos que dos series se relaciona o que puede haber alguna relación entre ellas, así que al no saber si se relaciona con la serie con retardo se puede hacer un scatter plot para explorar esta posible relación.

#Scatterplot de dos series y sus retrasos
rec<-scan("recruit.dat")
rec=ts(rec)
par(mfrow=c(3,3),mar=c(2.5,4,4,1))
for(h in 0:8){
      plot(lag(soi,-h),rec,main=paste("Soi(t-",h,")"),ylab="Rec(t)",xlab="",col="2")}

Scatterplot-SOI-RecAl cierre de la entrada ¿qué es una serie de tiempo?, se comentó que las serie SOI y Recruitment muestran una posible relación, mientras SOI se refiere al promedio de la temperatura en el océano pacífico, Recruitment cuenta con los datos de la estimación de nueva población de peces en el océano pacífico. Así que la gráfica anterior hace una exploración de la posible relación entre los datos con ciertos retrasos.

Otra herramienta para explorar la serie es analizar sus espectro de frecuencias, esto está relacionado con la versión de análisis de señales. Para lo cual se hace uso de la transformada de Fourier y se busca analizar las frecuencias  que dominan el comportamiento de la serie.

#Ejemplo de análisis espectral
t=1:500
x=2*cos(2*pi*t/50+.6*pi)
I=abs(fft(x)/sqrt(500))^2
P=(4/500)*I
f=0:250/500
plot(f,P[1:251],type="l",xlab="Frecuencias",main="Ejemplo de Espectro de Frencuencias",col="4")
abline(a=seq(0,.5,.02),lty="dotted")

 Espectro_de_Frecuencias

En la gráfica anterior se muestra la frecuencia dominante de la serie, esto permite empezar hacer un análisis de la señal. Pese a que los datos fueron generados en R project, uno puede recuperar información de la serie original al hacer el análisis espectral de la serie.

Algo de modelos de componentes no observables.

Antes de explicar en breve que es eso de los modelos de componentes no observables hago la descomposición de la serie SOI.

#Descomposición de la serie
plot(stl(ts(soi,freq=12),s.window="periodic"),col="4",main="Descomposición de SOI")

Descomposición_de_SOI

La gráfica muestra la serie, la parte estacional , la tendencia y los residuos. Cada una parece dar información sobre la serie. Lo que hacen los modelos de componentes no observables es tomar la idea de que cada señal está compuesta por cuatro componentes no observados:Tendencia, Estacionalidad, Ciclo e Irregularidad.

La idea  modelo es la siguiente:

Señal=f(T,S,C,I)

Donde f() es una función y tienen esas cuatro variables. Los modelos más sencillos y más usados son los que se llaman aditivos y multiplicativos. Con esto se indica que:

Señal=T+S+C+I

o

Señal=TxSxCxI

En series de tiempo no tan grandes la parte Cíclica es notoria, pero cuando no es suficientemente grande la serie la parte cíclica y la tendencia parecen altamente relacionadas y se opta por considerar solo un componente Tendencia-Ciclo. Así se reduce el modelo a:

Señal=T+S+I

o

Señal=TxSxI

Por lo cual en la gráfica anterior se mostraban solo 3 componentes. Entonces lo que resta es analizar como extraer información de la tendencia y de la estacionalidad.

La forma de elegir qué modelo usar es observando el comportamiento de la gráfica de la serie, si la estacionalidad se muestra contante con respecto  la tendencia el modelo es Aditivo, si la estacionalidad incremente con respecto a la tendencia el modelo es Multiplicativo.

#Serie aditiva
x=20*cos(2*pi*t/50+.6*pi)+rnorm(500,0,10)
Aditiva=x+t
plot.ts(Aditiva,col="2",main="Serie Aditiva",xlab="Tiempo",ylab="Aditiva")

#Ejemplo Multiplicativo

souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
plot.ts(souvenir,col="4",main="Ventas de tiendas de Recuerdos",xlab="Tiempo",ylab="Ventas")

 Series_Aditiva

Ejemplo_Multiplicativa

Los primeros datos son simulaciones de una serie aditiva y los otros datos son tomados desde la página de Rob Hyndman. Observamos en las gráficas como el comportamiento es distinto, mientras la primera mantiene la amplitud constante, la segunda se altera con respecto a la tendencia.

Lo que suele estudiar en textos al tener la noción de modelos de componentes no observables, se estudia como extraer la tendencia y como extraer la estacionalidad. Para cada una se usan ciertas técnicas que en otra entrada comentaré.

La parte de Irregularidad o ruido, es analizada con pruebas estadísticas para validar que cumplen ciertas condiciones que son convenientes e importantes teóricamente. Para los modelos ARIMA se requiere explicar y ejemplificar ciertos conceptos que están íntimamente ligados a la parte Irregular de la serie y que son lo que da origen a los modelos.

Espero que la entrada de una idea general de como explorar la serie.  Otro problema que surge es como analizar una familia de series de tiempo, para ellos se pueden usar técnicas como Análisis de Compones Principales (PCA) o Análisis de aglomerados (cluster).  Espero en otra entrada hablar al respecto.

Referencias:

1.-http://anson.ucdavis.edu/~shumway/

2.-https://www.otexts.org/fpp/

3.-http://www.statsoft.com/Textbook/Time-Series-Analysis

4.-http://www.medellin.unal.edu.co/~ndgirald/members.htm

5.-http://pandas.pydata.org/pandas-docs/stable/timeseries.html

6.-Analysis of financial time series

7.-Nonlinear time series analysis

8.-Econometric Models Economic Forecasts

¿Qué es una serie de tiempo?

Antes de intentar contestar la pregunta de la entrada, comparto 4 gráficos.
Temperatura_MediaSpeech_recordingNYSE_returnsCaminata_Aleatória

 

Las 4 gráficas anteriores muestras series de tiempo, de un modo muy informal las series de tiempo son datos con la etiqueta del tiempo. Ejemplo, el valor de la temperatura de la ciudad con el registro de la hora en la cual fue medida o el valor de la tasa de cambio registrada junto con el dato del día en que se observó ese valor. Claro con la idea de que la medición del tiempo es regular, se realiza cada hora o cada día o cada semana.

Las series de tiempo son una herramienta fácil para hacer un modelo y pronósticos a corto plazo. La ventaja es que hay muchos datos en la red que pueden ser estudiados como series de tiempo, que van desde indicadores económicos, sísmicos, datos financieros, etc.

La cantidad de técnicas, libros, artículos y campos de estudio de las series de tiempo es amplio, basta con echarle un vistazo en la red y sobran fuentes de donde aprender sobre el tema. La perspectiva también es distinta, por un lado desde el lado de matemáticas suelen ser aspectos teóricos, desde física suelen ser métodos no siempre estándar y la perspectiva de finanzas o economía suelen ser métodos que suelen tener mucha información y versan entre teoría y ejemplos prácticos, esto último debido a que muchos datos económicos son estudiados como series de tiempo.

El motivo de escribir estas entradas es compartir ejemplos y la metodología para hacer modelos lineales, o ARIMA y un poco sobre análisis espectral de una serie. Pero principalmente dar ejemplos.

Todas las explicaciones no serán exhaustivos, pero poco a poco les dejo texto a referencias para que se consulte con mayor detalle.

El software que uso en las entrada es R project, pero bien se puede usar o replicar los ejemplos con Gretl o con Eviews .

En Python se tienen las librerías de Pandas y en C++ hay varias librerías para manejo de series de tiempo, si les es más como usar estos dos lenguajes de programación está perfecto.

La única ventaja que veo con R es que si uno no sabe programar, basta copiar y pegar el código en su computadora para replicar los ejemplos y posteriormente podrá migrar a otro software más robusto con mayor facilidad, como Python.

Los datos de los ejemplos se pueden descargar de las páginas de Robert S. Shumway y de Ruey S.Tsay. Las fuentes de algunos datos las iré comentando en cada entrada.

Series de tiempo

Las series de tiempo son una herramienta fácil para hacer pronósticos a corto plazo cuando no se tiene mucha información sobre el sistema o el problema que se analiza.

Los principales modelos de series de tiempo son Medias Móviles (Moving Averages), Autorregresivos (Autoregressions)  y combinación de ellos  Autorregresivos con Medias Móviles (ARMA). Las herramientas estadísticas básicas que se usan son la autocorrelación y la autocorrelación parcial. Estas son las herramientas básicas para hacer algunos pronósticos no tan malos.

Pero la familia de modelos y técnicas para hacer predicciones por medio de series de tiempo van desde técnicas de Maquina de Soporte Vectorial, Redes Neuronales,  Cadenas de Markov Ocultas, Técnicas Lineales (Filtro de Kalman).

Faltaría agregas a todas las técnicas anteriores los métodos no lineales y las técnicas desarrolladas desde sistemas complejos para analizar series de tiempo, que son cercanas a las técnicas desarrolladas por los geofísicos.

Pero en las entradas solo hablaré de los métodos básicos y sobre todo de ejemplos de como analizar las series. Cada una de las técnicas anteriores requieren mucho más que unas cuantas entradas para explicarlas o hablar de ellas.

Hacemos unos ejemplos con datos generados en R para mostrar lo que es un ruido blanco, una caminata aleatoria, una serie generada por un proceso de medias móviles y procesos autoregressivo.

#Ruido blanco
ruido=rnorm(1000,0,1)
plot.ts(ruido)

#Otro modo de hacer la gráfica es definiendo la series de tiempo
ruido<-ts(ruido)
plot(ruido,col="2")

#Simulación de Series Media Movil

ma<-filter(ruido,rep=c(1,6)/6,sides=2)
plot.ts(ma,main="Media Móvile",col="4")

#Simulación de autoregresivo

rg<-filter(ruido,filter=c(1,-.9),method="recursive")
plot.ts(rg,col="4",main="Autoregresivo")

#Caminata Aleatoria

ruidoW=ruido+.2;x=cumsum(ruido);xd=cumsum(ruidoW)
plot.ts(xd,col="2",main="Caminata Aleatoria")
lines(x,col="5")
lines(.2*(1:1000),lty="dashed",col="3")

Ruido_BlancoMAAutoregresivoCaminata_Ale

Hacemos el cálculo de la autocorrelación y de la autocorrelación parcial para las dos series simuladas, es decir, el modelo de Medias Móviles y el Autorregresivo.

#Cálculo de la autocorrelaciones
acf(ma[3:997],main="Modelo de Medias Móviles")
acf(rg,main="Modelo de Autoregresivo")

#Cálculamos las autocorrelaciones parciales

pacf(ma[3:997],main="Modelo de Medias Móviles")
pacf(rg,main="Modelo de Autorregresivo")

 ACF_MAPacf_MAACF_AutorregresivoPacf_Autorregresivo

Estas últimas gráficas dan información de como ir determinando un modelo ARMA, donde se determina la parte generada por un proceso de medias móviles y otra parte por un autorregresivo. Con un poco de cuidado se observa que el número de barras negras de la primera gráfica corresponde con el número de elementos considerados para la media móvil y el número de barras negras en la última gráfica corresponde con el número de elementos para la regresión del modelo autorregresivo.

Es por ellos la importancia de estas dos funciones, que se pueden calcular en R project prácticamente con una función base.

Hago la gráfica de datos reales y se observará que tiene similitud con el modelo de Media Móvil con la única diferencia de que los datos tendrán una tendencia. Estos datos corresponde al promedio de nacimientos en la ciudad de New York desde Enero de 1946 hasta Diciembre de 1956.

#Datos del promedio de nacimientos en NY

births <- scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
seriedetiempo_births<- ts(births, frequency=12, start=c(1946,1))

plot(seriedetiempo_births,main="Porcentaje de nacimientos", xlab="Años",ylab="Porcentaje de nacimientos",col="2")

 Porcentaje_nacimientos_NY

Lo que trato de mostrar con esta última gráfica es que los modelos ARMA son un buen modelo, para analizar cierto tipo de datos. Además en el código se observa que no se descargan los datos , sino se toman desde la página de Rob J Hyndman, esto se puede hacer desde la función scan y todas las funciones que la requieren, como read.table.

Pero la parte de “analizar datos” o las series, es análogo a pensar que tenemos una “señal ruidosa”, es decir; busca desarrollar métodos para extraer la señal y eliminar el ruido que tienen. Esa es la idea en general del análisis de series de tiempo.

Cómo ejemplo del párrafo anterior hago un pequeño  ejemplo.

#Señal y ruido
t=1:750
c=2*cos(2*pi*t/50+.6*pi)
ruido2=rnomr(750)
par(mfrow=c(3,1))
plot.ts(c, main="Ejemplo de Señal y Ruido")
plot.ts(c+ruido2,col="2")
plot.ts(c+5*ruido2,col="4")

 Señal_y_ruido

La gráfica muestra una onda y como al agregarle un ruido blanco empieza a tener una forma parecida a datos que aparecen de la naturaleza o en cierto tipo de fenómenos.

Pensar en analizar una sola serie de tiempo, esta bien; pero puede presentarse o uno se puede preguntar si dos series tiene alguna relación y más aún si los fenómenos que las generan pueden parecer relacionados. Como ejemplo, es que exista una posible relación de la temperatura de la superficie del océano Pacífico Central con datos de la población de peces,  se revisan sus gráficas.

#Datos de SOI y de Poblaciones de Peces
#Considerando que ya se cargaron los datos en R o que se está en el directorio de los tatos

soi<-ts(soi,start=c(1950,1),frequency=12)
newfish<-ts(recruiment,star=c(1950,1),frequency=12)

par(mfrow=c(2,1))
plot(soi,main="Temperatura del Pacífico Central")
plot(newfish,main="Comportamiento de la Población de Peces")

 Dos_series_SOI-NF

La gráfica muestra que no es descabellado pensar que exita una relación entre las series de tiempo. Entonces para analizar su comportamiento lo que se hacemos uso de ACF y de la Correlación -Cruzada o Covarianza de las series.

#Relación entre series
par(mfrow=c(3,1))
acf(soi,50)
acf(newfish,50)
ccf(soi,newfish,50)

 Relación_entre_series

En otras entradas explico más respecto a lo que se hace con el manejo de varias series. En esta entrada solo pretendía mostrar algunos ejemplos y hablar un poco de lo que se requiere para hacer el análisis de las series y como elegir un modelo ARIMA para describir la serie y hacer algunos pronósticos.

 Referencias:

1.-http://anson.ucdavis.edu/~shumway/

2.-https://www.otexts.org/fpp/

3.-http://www.statsoft.com/Textbook/Time-Series-Analysis

4.-http://www.medellin.unal.edu.co/~ndgirald/members.htm

5.-http://pandas.pydata.org/pandas-docs/stable/timeseries.html

6.-Analysis of financial time series

7.-Nonlinear time series analysis