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
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")
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")
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")
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")
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")
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")
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
dlegorreta 
Debe estar conectado para enviar un comentario.