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

Anuncios

Un comentario sobre “Series con tendencia y estacionalidad.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s