Series de Tiempo

Introduccion al Analisis de Series de Tiempo

Librerias para series de tiempo

library(lubridate)

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(tseries)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

    ‘tseries’ version: 0.10-49

    ‘tseries’ is a package for time series analysis and
    computational finance.

    See ‘library(help="tseries")’ for details.

library(forecast)
library(ggplot2)
library(zoo)

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(chron)

Attaching package: ‘chron’

The following object is masked from ‘package:tseries’:

    is.weekend

The following objects are masked from ‘package:lubridate’:

    days, hours, minutes, seconds, years

library(dygraphs)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(readxl)
library(highcharter)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Fechas

Dependiendo el valor de entrada tenemos que usar alguna de las tres, el valor de salida siempre sera: YYYY-MM-DD

ymd(19931123)

[1] "1993-11-23"

dmy(23111993)

[1] "1993-11-23"

mdy(11231993)

[1] "1993-11-23"

Fechas y tiempo

mytimepoint = ymd_hm('1993-11-23 11:23', tz = "America/Argentina/Buenos_Aires")

class(mytimepoint)

[1] "POSIXct" "POSIXt" 

Para extraer los componentes:

minute(mytimepoint)

[1] 23

day(mytimepoint)

[1] 23

hour(mytimepoint)

[1] 11

year(mytimepoint)

[1] 1993

month(mytimepoint)

[1] 11

wday(mytimepoint, label = T, abbr = F)

[1] martes
Levels: domingo < lunes < martes < miércoles < jueves < viernes < sábado

Calculando intervalo entre dos momentos temporales:

time1 = ymd_hm('1993-11-23 11:23', tz = "America/Argentina/Buenos_Aires")


time2 = ymd_hm('1995-12-23 11:23', tz = "America/Argentina/Buenos_Aires")


myinterval = interval(time1, time2)


class(myinterval)

[1] "Interval"
attr(,"package")
[1] "lubridate"

Outliers y NA

mydata = read.csv('Rmissing.csv')

myts = ts(mydata$mydata)

plot(myts)

Podemos tratar los datos faltantes:

Nos llena los NA con la observación anterior

myts.NAlocf = na.locf(myts)

Rellena con el valor que queremos:

# myts.NAfil = fill(myts, 33)

Detección de outliers:

myts1 = tsoutliers(myts)
myts1

$index
[1]  45  98 172 211

$replacements
[1] 28.41242 39.78616 33.59838 37.96883

Limpia los NA y Outliers:

Se reemplazan los faltantes y se nivelan los outliers

# 
# mytsclean = tsclean(myts)
# 
# plot(mytsclean)

plot(tsclean(myts))

Datos de Variables Financieras:

Se importan datos de series temporales de los precios de Starbucks y Microsoft

Especificamos la fecha inicial, final y la frecuencia.

La seriede tiempo sbux inicia en:  1993 3 y finaliza en:  2008 3 con una frecuencia de:  12 

La seriede tiempo msft inicia en:  1993 3 y finaliza en:  2008 3 con una frecuencia de:  12

Subconjunto de una Serie Temporal:

subconjunto = window(sbux.ts, start = c(1993,3), end = c (1993,8))
subconjunto

      Mar  Apr  May  Jun  Jul  Aug
1993 1.13 1.15 1.43 1.46 1.41 1.44

Combinar Series de Tiempo:

combinadas_ts = cbind(sbux.ts, msft.ts)

Plot de Starbucks:

Plot de las dos Series:

Utilizando ZOO

Creamos una fecha con secuencias:

td = seq(as.Date('1993/3/1'),as.Date('2008/3/1'), 'months')

Otra alternativa es agarrar la fecha del data frame y convertirla:

td2 = as.Date(sbux$Date, format = '%m/%d/%Y')

Juntamos ambas Series:

# Combinamos el indice de tiempo a las dos series

sbux.z = zoo(x = sbux$Adj.Close, order.by = td)
msft.z = zoo(x = msft$Adj.Close, order.by = td)

# Extraer indices y precios

indices_sbux.z = index(sbux.z )


valores_sbux.z = coredata(sbux.z )


# Combinamos las dos series

combinadas_zoom = cbind(sbux.z, msft.z)

Importar datos directamente desde ZOO

# Es identico al creado anteriorimente

sbuxzoo = read.zoo('sbuxPrices.csv', 
                   format = '%m/%d/%Y',
                   sep = ",",
                   header = T)

Importar datos directamente de Yahoo Finance:

time series ends   2022-01-05
time series ends   2022-01-05

Grafico en conjunto:

Libreria dygraphs

Esta libreria contiene visualizaciones interactivas:

Ruido Blanco

Es un tipo especial de serie temporal, donde los datos no siguen ningún patrón.

Para que una serie sea ruido blanco tiene que tener una media constante, una varianza constante y no tener autocorrelacion en ningún periodo.

La autocorrelacion mide cuan correlacionada es una serie con versiones anteriores de si misma.

La falta de autocorrelacion es que no hay una relación clara entre valores pasados y valores presentes

\[p=corr(X_t,X_t-1) \]

Caminata Aleatoria

Es un tipo especial de serie de tiempo en donde los valores tienden a persistir en el tiempo y las diferencias entre periodos son simplemente ruido blanco.

Cualquier caminata aleatoria: Su precio es igual al precio anterior mas un ruido blanco

Las mejores estimaciones para los precios de hoy serán los precios de ayer

Estacionariedad

Prueba Dickey Fuller

La prueba nos sirve para testear si hay o no autocorrelacion en una serie temporal.

Esta propiedad hace referencia a cuando una serie es estable a lo largo del tiempo, es decir, cuando la media y la varianza son constantes y ademas no presenta tendencias

Un ejemplo de una serie débilmente estacionaria es una serie de ruido blanco, ya que la autocorrelacion entre dos puntos siempre es cero

Serie Estacionaria

Simulamos datos normalmente distribuidos de 1000 observaciones y aleatorios, el dickey fuller test lo vamos a ejecutar para ver si hay significancia o no. El p valor es menor a un nivel de significancia del 0.05, por lo tanto, la hipotesis nula es rechazada, por lo tanto la serie sera estacionaria.

    Augmented Dickey-Fuller Test

data:  x
Dickey-Fuller = -11.557, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary

Serie No Estacionaria

    Augmented Dickey-Fuller Test

data:  y
Dickey-Fuller = -2.6997, Lag order = 9, p-value = 0.2821
alternative hypothesis: stationary

Descomponer Serie

Warning in adf.test(nottem) : p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  nottem
Dickey-Fuller = -12.998, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary

Autocorrelacion

La correlación nos dice la similitud con respecto a las series temporales cambian sus valores. Si estamos con una sola serie temporal la correlación se calcula con valores pasados de la misma serie, es decir, la correlación entre una secuencia y si misma. Mide el nivel de semejanza entre una secuencia de varios periodos atrás y los datos reales, esta secuencia del pasado se llama retraso o rezago.

La función de autocorrelacion ACP Function nos va a proporcionar la autocorrelacion para cualquier retraso que consideremos.

Autocorrelacion Parcial

Son efectos de segunda mano, donde los valores presentes son afectados por valores pasados y estos mismos también son afectados por valores pasados

Autocorrelacion para ruido blanco

Autocorrelacion para ruido blanco.(Serie Estacionaria)

acf(x, plot = T)

Estacionalidad

Sugiere que ciertas tendencias aparecerán de forma cíclica, por ejemplo, las temperaturas suben y bajan a lo largo del tiempo. La series se pueden descomponer bajo efecto tendencia donde podemos quitar un patrón sostenido en los datos, la estacional los efectos cíclicos y el residual el error entre los datos reales y el modelo que estamos ejecutando.

La descomposición clásica puede ser aditiva o multiplicativa, en la aditiva podemos suponer que en cualquier momento del tiempo el valor observado es la suma de la tendencia, el objeto estacional y el efecto residual para ese periodo.

La descomposición multiplicativa nos dice que la serie original es un producto de los tres efectos anteriories.

Suavizado con SMA

Utilizaremos un Data Set y crearemos una variable de SMA de orden 3 y de orden 9.

Método Exponencial

etsmodel = ets(nottem)
plot(nottem, lwd = 3)
lines(etsmodel$fitted, col = "red")

# Forecast de 12 meses
# Intervalo de predicción: nivel de confianza 95%
plot(forecast(etsmodel, h = 12, level = 95))

# Holt Winters multiplicativo
etsmodmult = ets(nottem, model ="MMM") # MMM todo multiplicativo #error, tendencia, estacionalidad. Defecto ZZZ

# Comparacion serie original con el fitted luego de aplicar el modelo multiplicativo
plot(nottem, lwd = 3)
lines(etsmodmult$fitted, col = "red")

etsmodmultMAM = ets(nottem, model ="MAM") 
plot(nottem, lwd = 3)
lines(etsmodmultMAM$fitted, col = "red")

compare_models = as.data.frame(nottem)

compare_models$fitted_MMM = etsmodmult$fitted

compare_models$fitted_MAM = etsmodmultMAM$fitted

# Error cuadrático medio

cat('El error cuadratico medio para tendencia multiplicativa es: ', sum((compare_models$x - compare_models$fitted_MMM) / nrow(compare_models)), 'y para tendencia aditiva: ', sum((compare_models$x - compare_models$fitted_MAM) / nrow(compare_models)))

El error cuadratico medio para tendencia multiplicativa es:  0.163262 y para tendencia aditiva:  0.02354155

Modelos Predictivos

¿Como seleccionar el mejor modelo?

1) Escoger el modelo

2) Dividir los datos de entrenamiento y de prueba

3) Ajustar el modelo a los datos de entrenamiento

4) Evaluar el modelo en los datos de prueba

5) Re-Ajustar el modelo con todos los datos

6) Predecir los datos futuros

Criterios de selección del modelo:

AIC = AKaike: Evalúa la colección de modelos que tenemos y va a estimar la calidad de cada modelo en relación con los modelos restantes, el criterio va a penalizar en base al numero de parámetros que tiene, es decir, en base a su complejidad.

BIC = Baysesian: Tiene criterios Bayesianos. Utilizar los dos criterios es lo recomendado.

Conceptos Claves:

Overfitting: Sobreajuste del modelo. Puede no funcionar de forma correcta ya que esta sobre ajustado el modelo

Los residuos del modelo deben parecerse a ruido blanco, es decir, estacionarios

Modelo Autorregresivo

La autocorrelacion nos va a permitir desarrollar este modelo que se basa en tener en cuenta valores pasados. Es un modelo lineal donde los valores del periodo actual son la suma de resultados pasados multiplicados por un coeficiente mas un error.

El coeficiente siempre tiene que estar entre - 1 y 1. No puede ser superior en valor absoluto

El residuo van a ser diferencias impredecibles, si hay un patrón se va a identificar en las variables del modelo, el residuo es la diferencia entre el valor real y el estimado

Modelo AR(1)

$$ \[\begin{aligned} Y_i &= c + \phi_i Y_{i-1} + \epsilon_i \end{aligned}\]

$$

Modelo AR(2)

$$ \[\begin{aligned} Y_i &= c + \phi_i Y_{i-1} + \phi_i Y_{i-2} + \epsilon_i \end{aligned}\]

$$ Modelo AR(n)

$$ \[\begin{aligned} Y_i &= c + \phi_i Y_{i-1} + \phi_2 Y_{i-2} + ..\phi_n Y_{i-n} + \epsilon_i \end{aligned}\]

$$ Estos modelos de muchos procesos AR pueden contener ruido. Es importante primero verificar de forma teórica para no hacerlo mas complejo.

Las funciones de autocorrelacion y autocorrelacion parcial nos van a decir cuantos rezagos podemos incluir en el modelo.

ARMA

Estos modelos combinan los anteriores vistos. Toma en cuenta los valores pasados y los errores.

$$ \[\begin{aligned} Y_i &= c + \phi_i Y_{i-1} + \alpha_i \epsilon_{i-1} + \epsilon_i \end{aligned}\]

$$ Donde el coeficiente AR es:

$$ \[\begin{aligned} \phi_i \end{aligned}\]

$$ Y el coeficiente MA es:

\[ \begin{aligned} \alpha_i \end{aligned} \]

ARIMA

Cuando se trata de series no estacionarias utilizaremos el modelo autorregresivo integrado de medias moviles. La parte integrada del modelo explica el numero de diferencias no estacionales que debemos examinar para establecer estacionariedad.

ARIMA (p,d,q)

El modelo va a tener tres ordenes, p,d y q. Donde p = Componente AR, d = Orden de Integración y 1 = Componente MA.

Donde un ARIMA(p,0,q) = ARMA(p,q).

Numero de linces capturados anualmente en EEUU

En este gráfico obtenemos las funciones de autocorrelacion y el mismo gráfico anterior.

Observamos cierto comportamiento cíclico, no parece haber estacionalidad.

En la FACP podemos ver que el primer y segundo retraso son significativos.

Vemos informacion del data set:

tail(lynx)

Time Series:
Start = 1929 
End = 1934 
Frequency = 1 
[1]  485  662 1000 1590 2657 3396

Probamos un modelo AR(2)

modelarima = arima(lynx, order = c(2,0,0))

modelarima

Call:
arima(x = lynx, order = c(2, 0, 0))

Coefficients:
         ar1      ar2  intercept
      1.1474  -0.5997  1545.4458
s.e.  0.0742   0.0740   181.6736

sigma^2 estimated as 768159:  log likelihood = -935.02,  aic = 1878.03

Observamos los residuos del modelo:

Ejecutamos un MA(2):

Call:
arima(x = lynx, order = c(0, 0, 2))

Coefficients:
         ma1     ma2  intercept
      1.1407  0.4697  1545.3670
s.e.  0.0776  0.0721   224.5215

sigma^2 estimated as 855092:  log likelihood = -941.03,  aic = 1890.06

Hacemos un test de p valor para ver si la serie es estacionaria al 5 % de nivel de significancia:

La serie es estacionaria con un p valor de:  0.01

Probamos un modelo AR(4) ya que en la FACP se puede ver un rezago 4 significativo:

Podemos ver residuos dentro de la zona de no significación, que siguen aproximadamente una distribución normal y la serie temporal de los residuos vemos que es estacionaria para el modelo de orden 4.

Series: lynx 
ARIMA(4,0,0) with non-zero mean 

Coefficients:
         ar1      ar2     ar3      ar4       mean
      1.1246  -0.7174  0.2634  -0.2543  1547.3859
s.e.  0.0903   0.1367  0.1361   0.0897   136.8501

sigma^2 estimated as 748457:  log likelihood=-931.11
AIC=1874.22   AICc=1875.01   BIC=1890.64

    Ljung-Box test

data:  Residuals from ARIMA(4,0,0) with non-zero mean
Q* = 13.201, df = 5, p-value = 0.02157

Model df: 5.   Total lags used: 10

Haciendo Pronósticos:

El area oscura es el intervalo de confianza

Valores estimados:

arimaforecast$mean

Time Series:
Start = 1935 
End = 1949 
Frequency = 1 
 [1] 2980.7782 2114.6447 1361.7211  839.0137  668.7873  874.3079 1281.3753 1679.8363 1933.3503 1987.5494 1868.0533 1660.2155 1462.0024
[14] 1342.9317 1326.8693

Volatilidad

Hay modelos especiales que miden la volatilidad, es decir, la magnitud de los residuos, cuan volátiles sean las predicciones

Los modelos ARCH (Autorregresivos. con heterocedasticidad condicional) son los mas comunes para estudiar series temporales con volatilidad. A diferencia de la familia ARIMA este modelo tiene varias ecuaciones, una para la media y otra para la varianza.

La combinaron entre la heterocedasticidad y la condicional nos dice que la varianza dependerá de otros valores, es decir, de valores autorregresivos. Usaremos valores pasados para medir la varianza del periodo actual.

Predicciones

 ARIMA(2,0,2)(1,1,1)[12] with drift         : Inf
 ARIMA(0,0,0)(0,1,0)[12] with drift         : 33.18552
 ARIMA(1,0,0)(1,1,0)[12] with drift         : 21.24655
 ARIMA(0,0,1)(0,1,1)[12] with drift         : 4.650892
 ARIMA(0,0,0)(0,1,0)[12]                    : 31.10897
 ARIMA(0,0,1)(0,1,0)[12] with drift         : 30.61776
 ARIMA(0,0,1)(1,1,1)[12] with drift         : Inf
 ARIMA(0,0,1)(0,1,2)[12] with drift         : Inf
 ARIMA(0,0,1)(1,1,0)[12] with drift         : 22.19067
 ARIMA(0,0,1)(1,1,2)[12] with drift         : Inf
 ARIMA(0,0,0)(0,1,1)[12] with drift         : 3.449011
 ARIMA(0,0,0)(1,1,1)[12] with drift         : Inf
 ARIMA(0,0,0)(0,1,2)[12] with drift         : Inf
 ARIMA(0,0,0)(1,1,0)[12] with drift         : 21.62052
 ARIMA(0,0,0)(1,1,2)[12] with drift         : Inf
 ARIMA(1,0,0)(0,1,1)[12] with drift         : 4.220857
 ARIMA(1,0,1)(0,1,1)[12] with drift         : -0.587426
 ARIMA(1,0,1)(0,1,0)[12] with drift         : 23.21427
 ARIMA(1,0,1)(1,1,1)[12] with drift         : Inf
 ARIMA(1,0,1)(0,1,2)[12] with drift         : Inf
 ARIMA(1,0,1)(1,1,0)[12] with drift         : 15.62801
 ARIMA(1,0,1)(1,1,2)[12] with drift         : Inf
 ARIMA(2,0,1)(0,1,1)[12] with drift         : -1.704778
 ARIMA(2,0,1)(0,1,0)[12] with drift         : 23.66859
 ARIMA(2,0,1)(1,1,1)[12] with drift         : Inf
 ARIMA(2,0,1)(0,1,2)[12] with drift         : Inf
 ARIMA(2,0,1)(1,1,0)[12] with drift         : 14.07389
 ARIMA(2,0,1)(1,1,2)[12] with drift         : Inf
 ARIMA(2,0,0)(0,1,1)[12] with drift         : -0.02401758
 ARIMA(3,0,1)(0,1,1)[12] with drift         : 0.151561
 ARIMA(2,0,2)(0,1,1)[12] with drift         : -2.250542
 ARIMA(2,0,2)(0,1,0)[12] with drift         : Inf
 ARIMA(2,0,2)(0,1,2)[12] with drift         : Inf
 ARIMA(2,0,2)(1,1,0)[12] with drift         : 14.78926
 ARIMA(2,0,2)(1,1,2)[12] with drift         : Inf
 ARIMA(1,0,2)(0,1,1)[12] with drift         : -2.392319
 ARIMA(1,0,2)(0,1,0)[12] with drift         : 23.28728
 ARIMA(1,0,2)(1,1,1)[12] with drift         : Inf
 ARIMA(1,0,2)(0,1,2)[12] with drift         : Inf
 ARIMA(1,0,2)(1,1,0)[12] with drift         : 13.22757
 ARIMA(1,0,2)(1,1,2)[12] with drift         : Inf
 ARIMA(0,0,2)(0,1,1)[12] with drift         : 2.159611
 ARIMA(1,0,3)(0,1,1)[12] with drift         : -0.5346
 ARIMA(0,0,3)(0,1,1)[12] with drift         : 2.35328
 ARIMA(2,0,3)(0,1,1)[12] with drift         : -0.1783977
 ARIMA(1,0,2)(0,1,1)[12]                    : -4.59822
 ARIMA(1,0,2)(0,1,0)[12]                    : 21.08953
 ARIMA(1,0,2)(1,1,1)[12]                    : Inf
 ARIMA(1,0,2)(0,1,2)[12]                    : Inf
 ARIMA(1,0,2)(1,1,0)[12]                    : 11.00814
 ARIMA(1,0,2)(1,1,2)[12]                    : Inf
 ARIMA(0,0,2)(0,1,1)[12]                    : 0.002026895
 ARIMA(1,0,1)(0,1,1)[12]                    : -2.719768
 ARIMA(2,0,2)(0,1,1)[12]                    : -4.515401
 ARIMA(1,0,3)(0,1,1)[12]                    : -2.778658
 ARIMA(0,0,1)(0,1,1)[12]                    : 2.592817
 ARIMA(0,0,3)(0,1,1)[12]                    : 0.1663591
 ARIMA(2,0,1)(0,1,1)[12]                    : -3.911039
 ARIMA(2,0,3)(0,1,1)[12]                    : Inf

 Best model: ARIMA(1,0,2)(0,1,1)[12]                    

Estimaremos 3 años:

Prediccion de Linces:

La serie es estacionaria con un p valor de:  0.01

Descomposición de la serie:

Podemos observar, en la serie estacionaria, en la función de autocorrelacion parcial tenemos el primer y segundo rezago muy significativos, luego el cuarto también..

Modelo Auto ARIMA:

myarima = auto.arima(lynx, trace = T, 
           stepwise = F, 
           approximation = T)

 Fitting models using approximations to speed things up...

 ARIMA(0,0,0) with zero mean     : 2080.721
 ARIMA(0,0,0) with non-zero mean : 2006.724
 ARIMA(0,0,1) with zero mean     : 1971.633
 ARIMA(0,0,1) with non-zero mean : 1918.454
 ARIMA(0,0,2) with zero mean     : 1923.633
 ARIMA(0,0,2) with non-zero mean : 1890.07
 ARIMA(0,0,3) with zero mean     : 1911.329
 ARIMA(0,0,3) with non-zero mean : 1888.462
 ARIMA(0,0,4) with zero mean     : 1904.649
 ARIMA(0,0,4) with non-zero mean : 1888.861
 ARIMA(0,0,5) with zero mean     : 1906.728
 ARIMA(0,0,5) with non-zero mean : 1885.904
 ARIMA(1,0,0) with zero mean     : 1934.293
 ARIMA(1,0,0) with non-zero mean : 1926.793
 ARIMA(1,0,1) with zero mean     : 1902.375
 ARIMA(1,0,1) with non-zero mean : 1890.554
 ARIMA(1,0,2) with zero mean     : 1902.539
 ARIMA(1,0,2) with non-zero mean : 1888.371
 ARIMA(1,0,3) with zero mean     : 1904.57
 ARIMA(1,0,3) with non-zero mean : 1889.878
 ARIMA(1,0,4) with zero mean     : 1906.532
 ARIMA(1,0,4) with non-zero mean : 1888.773
 ARIMA(2,0,0) with zero mean     : 1906.818
 ARIMA(2,0,0) with non-zero mean : 1878.042
 ARIMA(2,0,1) with zero mean     : 1903.406
 ARIMA(2,0,1) with non-zero mean : 1879.57
 ARIMA(2,0,2) with zero mean     : 1905.591
 ARIMA(2,0,2) with non-zero mean : 1876.453
 ARIMA(2,0,3) with zero mean     : Inf
 ARIMA(2,0,3) with non-zero mean : Inf
 ARIMA(3,0,0) with zero mean     : 1904.652
 ARIMA(3,0,0) with non-zero mean : 1881.03
 ARIMA(3,0,1) with zero mean     : 1906.551
 ARIMA(3,0,1) with non-zero mean : Inf
 ARIMA(3,0,2) with zero mean     : Inf
 ARIMA(3,0,2) with non-zero mean : 1879.019
 ARIMA(4,0,0) with zero mean     : 1907.789
 ARIMA(4,0,0) with non-zero mean : 1876.139
 ARIMA(4,0,1) with zero mean     : Inf
 ARIMA(4,0,1) with non-zero mean : 1877.57
 ARIMA(5,0,0) with zero mean     : 1906.868
 ARIMA(5,0,0) with non-zero mean : 1878.465

 Now re-fitting the best model(s) without approximations...




 Best model: ARIMA(4,0,0) with non-zero mean 

Modelo Auto ARIMA con AR(8):

myarima = auto.arima(lynx, trace = T, 
           stepwise = F, 
           approximation = T,
           max.order=8, max.p=8)

 Fitting models using approximations to speed things up...

 ARIMA(0,0,0) with zero mean     : 2080.721
 ARIMA(0,0,0) with non-zero mean : 2006.724
 ARIMA(0,0,1) with zero mean     : 1971.633
 ARIMA(0,0,1) with non-zero mean : 1918.454
 ARIMA(0,0,2) with zero mean     : 1923.633
 ARIMA(0,0,2) with non-zero mean : 1890.07
 ARIMA(0,0,3) with zero mean     : 1911.329
 ARIMA(0,0,3) with non-zero mean : 1888.462
 ARIMA(0,0,4) with zero mean     : 1904.649
 ARIMA(0,0,4) with non-zero mean : 1888.861
 ARIMA(0,0,5) with zero mean     : 1906.728
 ARIMA(0,0,5) with non-zero mean : 1885.904
 ARIMA(1,0,0) with zero mean     : 1934.293
 ARIMA(1,0,0) with non-zero mean : 1926.793
 ARIMA(1,0,1) with zero mean     : 1902.375
 ARIMA(1,0,1) with non-zero mean : 1890.554
 ARIMA(1,0,2) with zero mean     : 1902.539
 ARIMA(1,0,2) with non-zero mean : 1888.371
 ARIMA(1,0,3) with zero mean     : 1904.57
 ARIMA(1,0,3) with non-zero mean : 1889.878
 ARIMA(1,0,4) with zero mean     : 1906.532
 ARIMA(1,0,4) with non-zero mean : 1888.773
 ARIMA(1,0,5) with zero mean     : Inf
 ARIMA(1,0,5) with non-zero mean : 1885.078
 ARIMA(2,0,0) with zero mean     : 1906.818
 ARIMA(2,0,0) with non-zero mean : 1878.042
 ARIMA(2,0,1) with zero mean     : 1903.406
 ARIMA(2,0,1) with non-zero mean : 1879.57
 ARIMA(2,0,2) with zero mean     : 1905.591
 ARIMA(2,0,2) with non-zero mean : 1876.453
 ARIMA(2,0,3) with zero mean     : Inf
 ARIMA(2,0,3) with non-zero mean : Inf
 ARIMA(2,0,4) with zero mean     : Inf
 ARIMA(2,0,4) with non-zero mean : Inf
 ARIMA(2,0,5) with zero mean     : Inf
 ARIMA(2,0,5) with non-zero mean : Inf
 ARIMA(3,0,0) with zero mean     : 1904.652
 ARIMA(3,0,0) with non-zero mean : 1881.03
 ARIMA(3,0,1) with zero mean     : 1906.551
 ARIMA(3,0,1) with non-zero mean : Inf
 ARIMA(3,0,2) with zero mean     : Inf
 ARIMA(3,0,2) with non-zero mean : 1879.019
 ARIMA(3,0,3) with zero mean     : Inf
 ARIMA(3,0,3) with non-zero mean : Inf
 ARIMA(3,0,4) with zero mean     : Inf
 ARIMA(3,0,4) with non-zero mean : Inf
 ARIMA(3,0,5) with zero mean     : Inf
 ARIMA(3,0,5) with non-zero mean : Inf
 ARIMA(4,0,0) with zero mean     : 1907.789
 ARIMA(4,0,0) with non-zero mean : 1876.139
 ARIMA(4,0,1) with zero mean     : Inf
 ARIMA(4,0,1) with non-zero mean : 1877.57
 ARIMA(4,0,2) with zero mean     : Inf
 ARIMA(4,0,2) with non-zero mean : Inf
 ARIMA(4,0,3) with zero mean     : Inf
 ARIMA(4,0,3) with non-zero mean : 1883.491
 ARIMA(4,0,4) with zero mean     : Inf
 ARIMA(4,0,4) with non-zero mean : Inf
 ARIMA(5,0,0) with zero mean     : 1906.868
 ARIMA(5,0,0) with non-zero mean : 1878.465
 ARIMA(5,0,1) with zero mean     : Inf
 ARIMA(5,0,1) with non-zero mean : Inf
 ARIMA(5,0,2) with zero mean     : Inf
 ARIMA(5,0,2) with non-zero mean : Inf
 ARIMA(5,0,3) with zero mean     : Inf
 ARIMA(5,0,3) with non-zero mean : Inf
 ARIMA(6,0,0) with zero mean     : 1904.575
 ARIMA(6,0,0) with non-zero mean : 1880.568
 ARIMA(6,0,1) with zero mean     : Inf
 ARIMA(6,0,1) with non-zero mean : Inf
 ARIMA(6,0,2) with zero mean     : Inf
 ARIMA(6,0,2) with non-zero mean : Inf
 ARIMA(7,0,0) with zero mean     : 1895.474
 ARIMA(7,0,0) with non-zero mean : 1881.193
 ARIMA(7,0,1) with zero mean     : 1878.376
 ARIMA(7,0,1) with non-zero mean : 1874.647
 ARIMA(8,0,0) with zero mean     : 1866.807
 ARIMA(8,0,0) with non-zero mean : 1860.933

 Now re-fitting the best model(s) without approximations...




 Best model: ARIMA(8,0,0) with non-zero mean 

Prediccion Arg

Prediccion saldos caja de ahorro ARG.

 ARIMA(2,1,2) with drift         : Inf
 ARIMA(0,1,0)            with drift         : 4048.257
 ARIMA(1,1,0) with drift         : Inf
 ARIMA(0,1,1)            with drift         : 4041.301
 ARIMA(0,1,0)                               : 4048.841
 ARIMA(1,1,1)            with drift         : 4033.881
 ARIMA(2,1,1) with drift         : Inf
 ARIMA(1,1,2)            with drift         : 4025.809
 ARIMA(0,1,2)            with drift         : 4023.928
 ARIMA(0,1,3)            with drift         : 4007.848
 ARIMA(1,1,3)            with drift         : 4006.809
 ARIMA(2,1,3) with drift         : Inf
 ARIMA(1,1,4) with drift         : Inf
 ARIMA(0,1,4) with drift         : Inf
 ARIMA(2,1,4) with drift         : Inf
 ARIMA(1,1,3)                               : 4019.882

 Best model: ARIMA(1,1,3)            with drift         

Filtro de Hampel

Para cada observación de la variable o serie de tiempo X, se calcula la mediana de una ventana (un subconjunto de valores). Ejemplo: Si nuestro conjunto de datos es una serie temporal de 100 observaciones diarias y nos fijamos en la observación 10, vamos a mirar un subconjunto de observaciones que se llamaran ventana, tenemos observaciones alrededor de donde estamos parados, vamos a tener vecinos, de cada lado de la observación tomada. 9 datos anteriores y 90 posteriores. Para cada ventana se calcula la mediana (medida representativa y robusta), es decir, si en esa ventana se encuentran los atípicos, la mediana va a ser una medida que represente el conjunto sin que este afectada por los atípicos. La media no es robusta porque se ve afectada por valores atípicos.

Luego se calculara para cada grupo de observaciones la distancia entre cada valor individual con la mediana del grupo, esto nos dará una desviación respecto a la mediana y para tener una medida representativa de todo el grupo se sacara la mediana de todas las desviaciones absolutas. Cambien conocida como mediana de la desviación absoluta con respecto a la mediana", que también se conoce como MAD (median absolute deviation).

La idea es que si una muestra difiere de la mediana en más de k desviaciones estándar, se considera un dato atípico y se reemplaza por el valor de la mediana.

Usualmente se suele seleccionar k=3, lo que se conoce como la regla de las “3 sigmas”. Pero esto puede depender del problema concreto.

Además entre el MAD y la desviación estándar se cumple la siguiente propiedad:

𝜎 ≈ 1.4826 MAD

Para el filtro de Hampel necesitamos definir dos cosas: 1. El tamaño de la ventana (cuántos vecinos vamos a considerar). 2. El número de desviaciones para identificar a los atípicos (k).

Un umbral más alto hace que el filtro sea más tolerante, uno más bajo identificará más puntos como valores atípicos.

library(pracma)

set.seed(33)


# Generar los datos aleatorios

x = numeric(1024)
z = rnorm(1024)
x[1] = z[1]
for (i in 2:1024) {
  x[i] = 0.4*x[i-1] + 0.8*x[i-1]*z[i-1] + z[i]
}


# Aplicar el Filtro de Hampel a nuestra serie temporal X con el parametro k que representa la ventana, es la cantidad de vecinos que vamos a considerar para hacer el calculo.

omad = hampel(x, k=20)


# Indice de los atípicos detectados:

omad$ind

 [1]  30  73  77  78  80  81  82  83 107 122 123 127 128 129 130
[16] 131 132 168 184 185 215 216 217 218 219 221 223 224 279 280
[31] 281 282 283 284 285 286 329 330 361 417 418 455 460 461 481
[46] 505 538 539 540 541 574 575 671 672 700 732 733 734 743 788
[61] 791 822 827 828 829 848 849 874 897 898 899 900 901 915 916
[76] 917 918 919 920 974

# Gráfico detectando atípicos
#
plot(1:1024, x, type="l")
points(omad$ind, x[omad$ind], pch=21, col="darkred")

# Nueva serie sin atípicos

x_new = omad$y


# Ambas series

plot(1:1024, x, type="l",col="red",xlim=c(0, 1050), ylim=c(-60,70))
par(new=TRUE)

plot(1:1024, x_new, type="l",col="blue",xlim=c(0, 1050), ylim=c(-60,70),axes= FALSE, xlab='', ylab='' )

NA
NA
NA