Introducción

Una serie de tiempo es una secuencia ordenada cronológicamente de observaciones y espaciadas entre sí de manera uniforme (usualmente los datos son independientes entre sí).

En la industria y el comercio, así como en la mayoría de las materias académicas cuantitativas, muchas variables se miden secuencialmente en el tiempo. Por ejemplo, en economía, las tasas de interés y las tasas de cambio se registran cada día. Se pueden medir diferentes variables en diferentes intervalos de tiempo, por ejemplo: la cantidad de lluvia capturada en un medidor se puede registrar para cada mes calendario; las temperaturas se pueden tomar cada hora en una estación meteorológica; o producto bruto nacional calculado anualmente. Cuando una variable se mide secuencialmente en el tiempo durante un intervalo fijo (o intervalo de muestreo), los datos forman una serie de tiempo.

El objetivo de del análisis de una serie de tiempo $x_t$ donde $t=1, \ldots, n$ es la realización de pronóstico.

Notación

A lo largo de este tema se representará a la serie de tiempo de tamaño $n$ con la notación: \[\left\{x_t : t=1,\ldots, n\right\}\], por lo que cada observación de la serie es representada como \[\left\{x_1,x_2, \ldots, x_n\right\}\].

La abreviación ${x_t}$ supondrá que se trabaja con una serie de tiempo de tamaño $n$ a menos de que se especifique. De igual forma que en regresión lineal $\hat{x}_t$ corresponde a los valores estimados o de predicción.

Algunos libros hacen notar que:

Si $t\le n$ entonces se le conocen a $\hat{x}_t$ como predicción (prediction)
Si $t > n$ entonces se le conocen a $\hat{x}_t$ como pronóstico (forecast)

Ejemplo: La siguiente base de datos hace referencia al consumo de energía eléctrica doméstica Total (Miles de millones de watts/hora) de manera mensual obtenida en el INEGI con fecha desde enero 2009 hasta diciembre 2017.

CONELDO=c(2861, 2833,   2686,   2798,   3123,   3232,   3377,   3545,   3558,   3275,   3079,   2792,   2825,   2769,   2655,   2735,   2996,   3282,   3594,   3544,   3633,   3422,   3220,   2859,   2998,   3062,   2907,   3012,   3349,   3537,   3688,   3632,   3792,   3608,   3274,   3049,   2999,   2990,   2983,   3068,   3279,   3440,   3636,   3670,   3835,   3605,   3436,   3126,   3126,   3120,   2905,   2992,   3272,   3524,   3732,   3807,   3793,   3631,   3481,   3203,   3230,   3176,   3005,   3119,   3382,   3543,   3820,   3985,   3950,   3743,   3512,   3211,   3229,   3212,   3043,   3188,   3483,   3678,   3896,   4142,   4190,   4063,   3847,   3482,   3415,   3330,   3063,   3284,   3677,   3964,   4211,   4296,   4310,   4142,   3858,   3591,   3371,   3331,   3297,   3396,   3753,   3967,   4322,   4341,   4398,   4213,   3867,   3509, 4000)

#Se convierte la base en un modelo de series de tiempo
CONELDO.ts=ts(CONELDO,start =2009,frequency = 12)
CONELDO.ts

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2009 2861 2833 2686 2798 3123 3232 3377 3545 3558 3275 3079 2792
## 2010 2825 2769 2655 2735 2996 3282 3594 3544 3633 3422 3220 2859
## 2011 2998 3062 2907 3012 3349 3537 3688 3632 3792 3608 3274 3049
## 2012 2999 2990 2983 3068 3279 3440 3636 3670 3835 3605 3436 3126
## 2013 3126 3120 2905 2992 3272 3524 3732 3807 3793 3631 3481 3203
## 2014 3230 3176 3005 3119 3382 3543 3820 3985 3950 3743 3512 3211
## 2015 3229 3212 3043 3188 3483 3678 3896 4142 4190 4063 3847 3482
## 2016 3415 3330 3063 3284 3677 3964 4211 4296 4310 4142 3858 3591
## 2017 3371 3331 3297 3396 3753 3967 4322 4341 4398 4213 3867 3509
## 2018 4000

class(CONELDO.ts)#El formato de la serie de tiempo es

## [1] "ts"

start(CONELDO.ts)#El Incio de la serie de tiempo es

## [1] 2009    1

frequency(CONELDO.ts)#La frecuencia de la serie de tiempo es

## [1] 12

summary(CONELDO.ts) #Información de la distribución de los datos

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2655    3120    3415    3443    3732    4398

Grafícanos la serie de tiempo para identificar patrones en el tiempo

plot.ts(CONELDO.ts, main='Consumo de energia domestica', xlab='Tiempo', ylab='Miles de millones de watts/hora' )

require(xts)

plot(as.xts(CONELDO.ts), major.ticks = 'quarters', main='Consumo de energia domestica')

En la gráfica de tiempo de los datos presenta características comunes a muchas series de tiempo. Por ejemplo, es evidente que el consumo de energía aumenta con el tiempo, por lo que parece haber una tendencia creciente. También hay un ciclo claro en los datos que tiene un período de un año, es decir, hay una clara variación estacional en los datos.

La comprensión de las causas probables de las características en la trama puede formar una guía útil al formular un modelo de serie de tiempo apropiado. En este caso, las posibles causas pueden incluir: aumento de la demanda de energía por el avance tecnológico; aumento de la demanda en ciertas épocas del año, durante períodos de vacaciones que podrían explicar la variación estacional. Alternativamente, una población en aumento podría explicar la tendencia cada vez mayor.

Un diagrama de serie temporal no solo es importante para revelar patrones y características de los datos, sino que también puede revelar valores atípicos potenciales o valores erróneos. Por ejemplo, los datos faltantes a veces se codifican con un valor negativo; obviamente, dichos valores tendrían que tratarse de manera diferente en el análisis, y ciertamente no querrían incluirse como observaciones al ajustar un modelo a los datos.

Para obtener una visión más clara de la tendencia, el efecto estacional se puede “eliminar” agregando los datos al nivel anual, lo que se puede lograr en R utilizando la función de aggregate. Se puede ver un resumen de los valores para cada temporada utilizando un diagrama de caja, con la función de cycle para extraer las estaciones para cada elemento de datos.

layout(1:2)
plot(aggregate(CONELDO.ts)); boxplot(CONELDO.ts ~ cycle(CONELDO.ts))

layout(1:1)

Ejemplo 2

La siguiente tabla contiene la información mensual del gasto en electricidad(en millones de KWh), la producción de cerveza y chocolate (en toneladas) en Australia en el periodo de enero de 1958 a diciembre de 1990, con información obtenida de the Australian Bureau of Statistics (ABS).

www = "D:/OneDrive - UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO/Documentos/ClasesFacultad/Curso 2021-2/Clase 04 - Introducción Series de Tiempo/cbe.dat.txt"
cbe = read.table(www, header=T)
head(cbe)

##   choc beer elec
## 1 1451 96.3 1497
## 2 2037 84.4 1463
## 3 2477 91.2 1648
## 4 2785 81.9 1595
## 5 2994 80.5 1777
## 6 2681 70.4 1824

Se convierte la información a un objeto de serie de tiempo, en el que se ingresará el mes y año inicial, así como mes-año final, además de ingresar la frecuencia de información

 elec.ts = ts(cbe[,3], start=1958, freq=12)
 beer.ts = ts(cbe[,2], start=1958, freq=12)
 choc.ts = ts(cbe[,1], start=1958, freq=12)

 plot(cbind(elec.ts, beer.ts, choc.ts),
main="Producción de chocolate, cerveza, y electricidad: 1958-1990")

Se muestra tendencias de estacionalidad y tendencia creciente en las tres variables, lo cual probablemente se deba al incremento de la población, así como sus exportaciones en Oceanía.

Se puede trabajar con las tres series al mismo tiempo, constituyendo una múltiple serie de tiempo. En R ayuda la función ts.intersect, lo que permite una serie de características importantes

#Se juntan las series 
 ap.elec = ts.intersect(elec.ts, beer.ts, choc.ts)
##Se obtiene fecha inicial
start(ap.elec);

## [1] 1958    1

##Se obtiene fecha final
end(ap.elec)

## [1] 1990   12

##Se tiene una sóla base datos
ap.elec[1:3, ]

##      elec.ts beer.ts choc.ts
## [1,]    1497    96.3    1451
## [2,]    1463    84.4    2037
## [3,]    1648    91.2    2477

#Se grafica de manera conunta
plot(ap.elec, main="Producción de chocolate, cerveza, y electricidad: 1958-1990")

##Se verifica tendencia
plot(aggregate(ap.elec))

Ejemplo pasajeros de avión

El siguiente ejemplo, es uno de más clásicos para el desarrollo de series de tiempo. La base se presenta en el texto de Box and Jenkins (1970), y es de gran utilidad para ilustrar la importancia de los conceptos de la exploración de los datos de series de tiempo.

La base contiene el número de personas (en miles) que viajan por mes en una aerolínea durante el periodo de 1949-1960.

data(AirPassengers)
AP = AirPassengers
AP

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432

plot(AP, main="Personas que viajan")

De la anterior gráfica se observa que Junio a Septiembre son los meses donde la gente viaja más lo que concuerda con los meses de periodo vacacional o de verano.

Se examina la tendencia y los valores que presentan observaciones atípicas, para su respectivo tratamiento.

layout(1:2)
plot(aggregate(AP)); boxplot(AP ~ cycle(AP))

Se observa que existe una tendencia a la alta en el número de vuelos conforme avanza el tiempo, lo que es una clara evidencia de que debe fortalecer esta industria ya que presenta una alza, además de que no hay demasiados valores atípicos por lo que implica una industria estable.

Correlación y Correlograma

Una herramienta útil para el desarrollo de series de tiempo es el correlograma, ya que permite identificar gráficamente las correlaciones.

Covarianza

La covarianza es el grado de variación conjunta entre una variable aleatoria $x$ y una variable aleatoria $y$ respecto a sus medias. Matemáticamente se define como: \[Cov(X,Y)=E\left[ \left( X-E[X]\right) \left( Y-E[Y]\right)\right],\] de esta manera se obtiene la asociación lineal entre las variables $(x,y)$: Si $x$ y $y$ tienden a incrementar juntas entonces se obtendrá una covarianza positiva, si $y$ decrece cuando $x$ crece entonces su covarianza será negativa. Si en las muestras aleatorias no se conocen su media se puede realizar un estimador de máxima verosimilitud para la covarianza dado por: \[Cov(x,y)=\frac{1}{n-1}\sum_1^n(x_i-\bar{x}_i) (x_{i-1}-\bar{y}_i),\] donde $\bar{x}_i=\frac{1}{n}\sum_1^n(x_i)$, $\bar{y}_i=\frac{1}{n}\sum_1^n(x_{i-1})$, cabe destacar, que algunos autores mencionan una simplificación de la anterior ecuación al sólo dividirlo entre $n$, sin embargo, esto es recomendable sólo cuando $n$ es un tamaño de muestra grande, pues así el efecto de $-1$ no aporta mucha significancia en el cálculo. ### Autocorrelación La correlación es la asociación lineal entre la variable $x$ y $y$ de una forma “estandarizada” a la covarianza. Esta es definida como: \[cor(X,Y)=\frac{cov(x,y)}{\sigma_x\sigma_y}\]

La correlación de Pearson tiene le ventaja de que toma valores entre -1 a 1 por lo que valores -1 implican una correlación negativa, es decir si la variable $x$ crece $y$ decrece en la misma proporción, mientras que 1 implica una correlación positiva perfecta, es decir si $x$ crece $y$ crece en prácticamente la misma proporción.

En series de tiempo se buscan observar si hay observaciones correlacionadas, dado que la serie sólo tiene una muestra aleatoria se realiza la autocorrelación con un retraso (lag) $k$ , así se calcula la autocorrelación: \[\rho_k=cor(x_{t}, x_{t+k})=\frac{cov(x_{t}, x_{t+k})}{\sqrt{Var^2(x_t)}},\]

La anterior cifra puede resumirse a: \[\rho_k=cor(x_{t}, x_{t+k})=\frac{\sum_{i=1}^{n-k} (x_i-\bar{x}_t)(x_{i+k}-\bar{x}_t)}{\sum_{i=1}^{n } (x_i-\bar{x}_t)^2},\] Note que: $\rho_0=1$. En R esto puede ser graficado con (acf(x))

acf(AP)

El anterior correlograma correponde a los pasajeros que vuelan de 1949-1960, el cual demuestra que las observaciones están altamente autocorrelacionadas entre sí.

Las líneas puntadas representan el 5% de significancia para probar que $\rho_k=0$

acf(CONELDO.ts)

El anterior correlograma correponde al consumo de energia domestica en México desde 2009 a 2017, el cual implica que ciertos valores son correlaccionados.

Estacionariedad y doble estacionaridad

Se dice que una serie de tiempo es estacionaria cuando es estable a lo largo del tiempo, es decir, la media y varianza son constantes en el tiempo, en ella los valores de la serie oscilan alrededor de la media $\mu_x=E[x_t]$ y la varianza $\gamma_k=cov(x_t, x_{tk})$, algunos autores llaman a la serie proceso de estacionariedad de segundo orden.

Descomposición de modelos

Modelo aditivo

El análisis de las series de tiempo se basa en la suposición de que los valores que toma la variable es consecuencia de 3 componentes cuya actuación conjunta da como resultados los valores medido

1. Componente de tendencia: Es el cambio de los valores a lo largo plazo que se produce con el paso de tiempo.
1. Componente estacional: Periodicidad o variación de la información que tiene un comportamiento similar en cierto periodo(mensual, semestral, anual, etc.)
1. Componente aleatorio: este componente no corresponde a ningún comportamiento, sino que es el resultado de factores aleatorios.

Así la serie es resultado de: \[x_t=m_t+s_t+\epsilon_t,\] donde. * $m_t$ componente tendencia al tiempo $t$. * $s_t$ componente estacional al tiempo $t$ * $\epsilon_t$ componente aleatorio al tiempo $t$

Este modelo es conocido como el modelo de descomposición clásico, este por defecto es el más óptimo para realizar series de tiempo, sin embargo hay otras tranformaciones ideales para cada tipo de datos, por ejemplo.

Si el efecto estacional tiende a aumentar a medida que aumenta la tendencia, un modelo multiplicativo puede ser más apropiado: \[x_t=m_t\times s_t+\epsilon_t,\] Si el modelo es múltiplicativo entonces se puede aplicar una transormación para para producir un módelo aditivo. \[log(x_t)=m_t+ s_t+\epsilon_t,\] Se requiere cierto cuidado cuando se aplica la transformación inversa exponencial, ya que el efecto suele ser sesgar las predicciones. Si las series residuales $ x_t $ se distribuyen normalmente con media 0 y varianza $ sigma ^ 2 $, entonces el valor predicho en el tiempo $t$ viene dado por \[\hat{x}_t=exp\{m_t+ s_t+\epsilon_t\},\]

Estamación de tendencia y estacionalidad

Hay varios métodos para estimar la tendencia $m_t$, sin embargo, lo más famosa es crear una serie de tiempo centrada en $x_t$. Lo más óptimo es realizar una media móvil de 12, 6 anteriores y 6 posteriores a $x_t$, de la forma:

\[\hat{m}_t=\frac{1/2 x_{t-6}+x_{t-5}+x_{t-4}+\ldots+x_{t}+\ldots+x_{t+5}+1/2x_{t+6}}{12}\]

Para estimar la estacionalidad ($S_t$) puede ser obtenido como: \[\hat{s}_t=x_t-\hat{m}_t\] Puede ser evidente, pero se supone que $\epsilon$ tiene media cero es por ello que no se considera en el cálculo. Po consiguiente en el modelo multiplicativo $\hat{s}_t=x_t/\hat{m}_t$

En R, esta estimación puede obtenerse a través del siguiente código:

#elect.ts es el consumo de energia en australia
plot(decompose(elec.ts, type = "additive"))

descomposicion=decompose(elec.ts, type = "additive")
descomposicion

## $x
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1958  1497  1463  1648  1595  1777  1824  1994  1835  1787  1699  1633  1645
## 1959  1597  1577  1709  1756  1936  2052  2105  2016  1914  1925  1824  1765
## 1960  1721  1752  1914  1857  2159  2195  2287  2276  2096  2055  2004  1924
## 1961  1851  1839  2019  1937  2270  2251  2382  2364  2129  2110  2072  1980
## 1962  1995  1932  2171  2162  2489  2424  2641  2630  2324  2412  2284  2186
## 1963  2184  2144  2379  2383  2717  2774  3051  2891  2613  2600  2493  2410
## 1964  2390  2463  2616  2734  2970  3125  3342  3207  2964  2919  2764  2732
## 1965  2622  2698  2950  2895  3200  3408  3679  3473  3154  3107  3052  2918
## 1966  2786  2739  3125  3033  3486  3661  3927  3851  3456  3390  3280  3166
## 1967  3080  3069  3340  3310  3798  3883  4191  4213  3766  3628  3520  3322
## 1968  3250  3287  3552  3440  4153  4265  4655  4492  4051  3967  3807  3639
## 1969  3647  3560  3929  3858  4485  4697  4977  4675  4596  4491  4127  4144
## 1970  4014  3994  4320  4400  5002  5091  5471  5193  4997  4737  4546  4498
## 1971  4350  4206  4743  4582  5191  5457  5891  5618  5158  5030  4800  4654
## 1972  4453  4440  4945  4788  5425  5706  6061  5846  5242  5408  5114  5042
## 1973  5008  4657  5359  5193  5891  5980  6390  6366  5756  5640  5429  5398
## 1974  5413  5141  5695  5554  6369  6592  7107  6917  6353  6205  5830  5646
## 1975  5379  5489  5824  5907  6482  6795  7028  6776  6274  6362  5940  5958
## 1976  5769  5887  6367  6165  6868  7201  7601  7581  7090  6841  6408  6435
## 1977  6176  6138  6717  6470  7312  7763  8171  7788  7311  6679  6704  6724
## 1978  6552  6427  7105  6869  7683  8082  8555  8386  7553  7398  7112  6886
## 1979  7077  6820  7426  7143  8261  8240  8977  8991  8026  7911  7510  7381
## 1980  7366  7414  7824  7524  8279  8707  9486  8973  8231  8206  7927  7999
## 1981  7834  7521  8284  7999  8940  9381 10078  9796  8471  8572  8150  8168
## 1982  8166  7903  8606  8071  9178  9873 10476  9296  8818  8697  8381  8293
## 1983  7942  8001  8744  8397  9115  9773 10358  9849  9083  9143  8800  8741
## 1984  8492  8795  9354  8796 10072 10174 11326 10744  9806  9740  9373  9244
## 1985  9407  8827  9880  9364 10580 10899 11687 11280 10208 10212  9725  9721
## 1986  9846  9407 10265  9970 10801 11246 12167 11578 10645 10613 10104 10348
## 1987 10263  9973 10803 10409 11458 11845 12559 12070 11221 11338 10761 11012
## 1988 10923 10790 11427 10788 11772 12104 12634 12772 11764 11956 11646 11750
## 1989 11485 11198 12265 11704 12419 13259 13945 13839 12387 12546 12038 11977
## 1990 12336 11793 12877 11923 13306 13988 14002 14338 12867 12761 12449 12658
## 
## $seasonal
##             Jan        Feb        Mar        Apr        May        Jun
## 1958 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1959 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1960 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1961 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1962 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1963 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1964 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1965 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1966 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1967 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1968 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1969 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1970 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1971 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1972 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1973 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1974 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1975 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1976 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1977 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1978 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1979 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1980 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1981 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1982 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1983 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1984 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1985 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1986 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1987 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1988 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1989 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
## 1990 -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
##             Jul        Aug        Sep        Oct        Nov        Dec
## 1958  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1959  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1960  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1961  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1962  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1963  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1964  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1965  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1966  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1967  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1968  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1969  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1970  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1971  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1972  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1973  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1974  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1975  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1976  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1977  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1978  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1979  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1980  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1981  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1982  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1983  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1984  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1985  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1986  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1987  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1988  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1989  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 1990  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 
## $trend
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1958        NA        NA        NA        NA        NA        NA  1703.917
## 1959  1772.958  1785.125  1797.958  1812.667  1830.042  1843.000  1853.167
## 1960  1936.500  1954.917  1973.333  1986.333  1999.250  2013.375  2025.417
## 1961  2071.375  2079.000  2084.042  2087.708  2092.833  2098.000  2106.333
## 1962  2194.958  2216.833  2236.042  2256.750  2278.167  2295.583  2312.042
## 1963  2438.583  2466.542  2489.458  2509.333  2525.875  2543.917  2561.833
## 1964  2708.458  2733.750  2761.542  2789.458  2814.042  2838.750  2861.833
## 1965  2989.125  3014.250  3033.250  3049.000  3068.833  3088.583  3103.167
## 1966  3194.750  3220.833  3249.167  3273.542  3294.833  3314.667  3337.250
## 1967  3473.500  3499.583  3527.583  3550.417  3570.333  3586.833  3600.417
## 1968  3734.917  3765.875  3789.375  3815.375  3841.458  3866.625  3896.375
## 1969  4079.000  4100.042  4130.375  4174.917  4210.083  4244.458  4280.792
## 1970  4506.500  4548.667  4586.958  4613.917  4641.625  4673.833  4702.583
## 1971  4848.417  4883.625  4908.042  4926.958  4949.750  4966.833  4977.625
## 1972  5082.750  5099.333  5112.333  5131.583  5160.417  5189.667  5228.958
## 1973  5413.792  5449.167  5492.250  5523.333  5546.125  5574.083  5605.792
## 1974  5841.792  5894.625  5942.458  5990.875  6031.125  6058.167  6067.083
## 1975  6157.875  6148.708  6139.542  6142.792  6153.917  6171.500  6200.750
## 1976  6406.792  6464.208  6531.750  6585.708  6625.167  6664.542  6701.375
## 1977  6901.417  6933.792  6951.625  6954.083  6959.667  6984.042  7011.750
## 1978  7190.583  7231.500  7266.500  7306.542  7353.500  7377.250  7405.875
## 1979  7589.000  7631.792  7676.708  7717.792  7755.750  7792.958  7825.625
## 1980  8013.708  8034.167  8041.958  8062.792  8092.458  8135.583  8180.833
## 1981  8423.083  8482.042  8526.333  8551.583  8576.125  8592.458  8613.333
## 1982  8769.250  8765.000  8758.625  8778.292  8793.125  8807.958  8803.833
## 1983  8822.833  8840.958  8875.042  8904.667  8940.708  8976.833  9018.417
## 1984  9345.083  9422.708  9490.125  9545.125  9593.875  9638.708  9697.792
## 1985  9947.542  9984.917 10024.000 10060.417 10094.750 10129.292 10167.458
## 1986 10384.000 10416.417 10447.042 10481.958 10514.458 10556.375 10599.875
## 1987 10866.833 10903.667 10948.167 11002.375 11059.958 11115.000 11170.167
## 1988 11400.208 11432.583 11484.458 11532.833 11595.458 11663.083 11717.250
## 1989 12125.625 12224.708 12295.125 12345.667 12386.583 12412.375 12457.292
## 1990 12748.625 12771.792 12812.583 12841.542 12867.625 12913.125        NA
##            Aug       Sep       Oct       Nov       Dec
## 1958  1712.833  1720.125  1729.375  1742.708  1758.833
## 1959  1865.625  1881.458  1894.208  1907.708  1922.958
## 1960  2034.458  2042.458  2050.167  2058.125  2065.083
## 1961  2116.208  2126.417  2142.125  2160.625  2176.958
## 1962  2328.750  2346.250  2364.125  2382.833  2406.917
## 1963  2583.708  2606.875  2631.375  2656.542  2681.708
## 1964  2881.292  2905.000  2925.625  2941.917  2963.292
## 1965  3111.708  3120.708  3133.750  3151.417  3173.875
## 1966  3363.250  3385.958  3406.458  3431.000  3453.250
## 1967  3616.583  3634.500  3648.750  3668.958  3699.667
## 1968  3924.292  3951.375  3984.500  4015.750  4047.583
## 1969  4314.167  4348.542  4387.417  4431.542  4469.500
## 1970  4725.417  4751.875  4777.083  4792.542  4815.667
## 1971  4991.667  5009.833  5026.833  5045.167  5065.292
## 1972  5261.125  5287.417  5321.542  5357.833  5388.667
## 1973  5642.833  5677.000  5706.042  5741.000  5786.417
## 1974  6080.167  6100.042  6120.125  6139.542  6152.708
## 1975  6233.583  6272.792  6306.167  6333.000  6366.000
## 1976  6728.792  6753.833  6781.125  6812.333  6854.250
## 1977  7039.458  7067.667  7100.458  7132.542  7161.292
## 1978  7444.125  7473.875  7498.667  7534.167  7564.833
## 1979  7862.417  7903.750  7936.208  7952.833  7973.042
## 1980  8204.792  8228.417  8267.375  8314.708  8370.333
## 1981  8643.083  8672.417  8688.833  8701.750  8732.167
## 1982  8798.583  8808.417  8827.750  8838.708  8831.917
## 1983  9074.417  9132.917  9174.958  9231.458  9288.042
## 1984  9737.250  9760.500  9806.083  9850.917  9902.292
## 1985 10209.917 10250.125 10291.417 10325.875 10349.542
## 1986 10640.833 10686.833 10727.542 10773.208 10825.542
## 1987 11231.708 11291.750 11333.542 11362.417 11386.292
## 1988 11757.667 11809.583 11882.667 11947.792 12022.875
## 1989 12517.542 12567.833 12602.458 12648.542 12715.875
## 1990        NA        NA        NA        NA        NA
## 
## $random
##               Jan          Feb          Mar          Apr          May
## 1958           NA           NA           NA           NA           NA
## 1959  322.8013238  431.7062717   61.1737196  345.8273655 -171.9981554
## 1960  283.2596571  436.9146050   90.7987196  273.1606988 -118.2064887
## 1961  278.3846571  399.8312717   85.0903863  251.7856988 -100.7898220
## 1962  298.8013238  354.9979384   85.0903863  307.7440321  -67.1231554
## 1963  244.1763238  317.2896050   39.6737196  276.1606988  -86.8314887
## 1964  180.3013238  369.0812717    4.5903863  347.0356988 -121.9981554
## 1965  131.6346571  323.5812717   66.8820530  248.4940321 -146.7898220
## 1966   90.0096571  157.9979384   25.9653863  161.9523655  -86.7898220
## 1967  105.2596571  209.2479384  -37.4512804  162.0773655  -50.2898220
## 1968   13.8429905  160.9562717  -87.2429470   27.1190321   33.5851780
## 1969   66.7596571   99.7896050  -51.2429470   85.5773655   -3.0398220
## 1970    6.2596571   85.1646050 -116.8262804  188.5773655   82.4185113
## 1971    0.3429905  -37.7937283  -14.9096137   57.5356988  -36.7064887
## 1972 -130.9903429  -19.5020616  -17.2012804   58.9106988  -13.3731554
## 1973   92.9679905 -152.3353950   16.8820530   72.1606988   66.9185113
## 1974   69.9679905 -113.7937283  -97.3262804  -34.3809679   59.9185113
## 1975 -280.1153429  -19.8770616 -165.4096137  166.7023655   50.1268446
## 1976 -139.0320095   62.6229384  -14.6179470  -18.2143012  -35.1231554
## 1977 -226.6570095 -155.9603950  -84.4929470  -81.5893012   74.3768446
## 1978 -139.8236762 -164.6687283  -11.3679470  -35.0476345   51.5435113
## 1979  -13.2403429 -171.9603950 -100.5762804 -172.2976345  227.2935113
## 1980 -148.9486762   19.6646050  -67.8262804 -136.2976345  -91.4148220
## 1981  -90.3236762 -321.2103950  -92.2012804 -150.0893012   85.9185113
## 1982 -104.4903429 -222.1687283   -2.4929470 -304.7976345  106.9185113
## 1983 -382.0736762 -200.1270616   19.0903863 -105.1726345 -103.6648220
## 1984 -354.3236762   12.1229384   14.0070530 -346.6309679  200.1685113
## 1985  -41.7820095 -518.0853950    6.1320530 -293.9226345  207.2935113
## 1986  -39.2403429 -369.5853950  -31.9096137 -109.4643012    8.5851780
## 1987 -105.0736762 -290.8353950    4.9653863 -190.8809679  120.0851780
## 1988   21.5513238   -2.7520616   92.6737196 -342.3393012 -101.4148220
## 1989 -141.8653429 -386.8770616  120.0070530 -239.1726345 -245.5398220
## 1990   86.1346571 -338.9603950  214.5487196 -516.0476345  160.4185113
##               Jun          Jul          Aug          Sep          Oct
## 1958           NA -630.1765408 -532.4981554   13.3442925   -2.4929470
## 1959 -316.2833116 -668.4265408 -504.2898220  -20.9890408   58.6737196
## 1960 -343.6583116 -658.6765408 -413.1231554    0.0109592   32.7153863
## 1961 -372.2833116 -644.5932075 -406.8731554  -50.9473741   -4.2429470
## 1962 -396.8666450 -591.3015408 -353.4148220  -75.7807075   75.7570530
## 1963 -295.1999783 -431.0932075 -347.3731554  -47.4057075   -3.4929470
## 1964 -239.0333116 -440.0932075 -328.9564887    5.4692925   21.2570530
## 1965 -205.8666450 -344.4265408 -293.3731554  -20.2390408    1.1320530
## 1966 -178.9499783 -330.5098741 -166.9148220   16.5109592   11.4237196
## 1967 -229.1166450 -329.6765408  -58.2481554   77.9692925    7.1320530
## 1968 -126.9083116 -161.6348741  -86.9564887   46.0942925   10.3820530
## 1969  -72.7416450 -224.0515408 -293.8314887  193.9276259  131.4653863
## 1970 -108.1166450 -151.8432075 -187.0814887  191.5942925  -12.2012804
## 1971  -35.1166450   -6.8848741  -28.3314887   94.6359592   31.0487196
## 1972   -8.9499783  -88.2182075  -69.7898220  -98.9473741  114.3403863
## 1973 -119.3666450 -136.0515408   68.5018446   25.4692925  -38.1596137
## 1974    8.5500217  119.6567925  182.1685113  199.4276259  112.7570530
## 1975   98.2166884  -93.0098741 -112.2481554  -52.3223741   83.7153863
## 1976   11.1750217  -20.6348741  197.5435113  282.6359592   87.7570530
## 1977  253.6750217  238.9901259   93.8768446  189.8026259 -393.5762804
## 1978  179.4666884  228.8651259  287.2101780   25.5942925  -72.7846137
## 1979  -78.2416450  231.1151259  473.9185113   68.7192925    2.6737196
## 1980   46.1333550  384.9067925  113.5435113  -50.9473741  -33.4929470
## 1981  263.2583550  544.4067925  498.2518446 -254.9473741  -88.9512804
## 1982  539.7583550  751.9067925 -157.2481554  -43.9473741 -102.8679470
## 1983  270.8833550  419.3234592  119.9185113 -103.4473741   -4.0762804
## 1984   10.0083550  707.9484592  352.0851780   -8.0307075  -38.2012804
## 1985  244.4250217  599.2817925  415.4185113  -95.6557075  -51.5346137
## 1986  164.3416884  646.8651259  282.5018446  -95.3640408  -86.6596137
## 1987  204.7166884  468.5734592  183.6268446 -124.2807075   32.3403863
## 1988  -84.3666450   -3.5098741  359.6685113  -99.1140408  101.2153863
## 1989  321.3416884  567.4484592  666.7935113 -234.3640408  -28.5762804
## 1990  549.5916884           NA           NA           NA           NA
##               Nov          Dec
## 1958  208.4771050  280.5773655
## 1959  234.4771050  236.4523655
## 1960  264.0604384  253.3273655
## 1961  229.5604384  197.4523655
## 1962  219.3521050  173.4940321
## 1963  154.6437717  122.7023655
## 1964  140.2687717  163.1190321
## 1965  218.7687717  138.5356988
## 1966  167.1854384  107.1606988
## 1967  169.2271050   16.7440321
## 1968  109.4354384  -14.1726345
## 1969   13.6437717   68.9106988
## 1970   71.6437717   76.7440321
## 1971   73.0187717  -16.8809679
## 1972   74.3521050   47.7440321
## 1973    6.1854384    5.9940321
## 1974    8.6437717 -112.2976345
## 1975  -74.8145616  -13.5893012
## 1976  -86.1478950  -24.8393012
## 1977 -110.3562283  -42.8809679
## 1978 -103.9812283 -284.4226345
## 1979 -124.6478950 -197.6309679
## 1980  -69.5228950   23.0773655
## 1981 -233.5645616 -169.7559679
## 1982 -139.5228950 -144.5059679
## 1983 -113.2728950 -152.6309679
## 1984 -159.7312283 -263.8809679
## 1985 -282.6895616 -234.1309679
## 1986 -351.0228950  -83.1309679
## 1987 -283.2312283   20.1190321
## 1988   16.3937717  121.5356988
## 1989 -292.3562283 -344.4643012
## 1990           NA           NA
## 
## $figure
##  [1] -498.75966 -639.83127 -150.13205 -402.49403  277.95649  525.28331
##  [7]  920.25987  654.66482   53.53071  -27.88205 -318.18544 -394.41070
## 
## $type
## [1] "additive"
## 
## attr(,"class")
## [1] "decomposed.ts"

De la anterior gráfica se observa la gráfica de valores observados, tendencia, estacionalidad, y aleatoriedad

#elect.ts es el consumo de energia en australia
plot(decompose(elec.ts, type = "mult"))

Vemos que las gráficas son similares, sin embargo, el cambio más notorio es que se presenta una mayor cantidad de error en esta última gráfica, pues el componente aleatorio presenta más alteraciones.

En R se puede seleccionar de manera individual la tendencia y estacionalidad de la siguiente forma

elec.decom = decompose(elec.ts, type='multiplicative')
trend = elec.decom$trend

seasonal = elec.decom$seasonal
 ts.plot(cbind(trend, trend * seasonal), lty=1:2)

 plot(elec.decom)

Pronósticos

El objetivo principal del pronóstico es predecir un valor futuro $x_{n+k}$ dada la historia pasada $\{x_1,x_2,\cdots,x_n\}$ de observaciones hasta el tiempo $n$. Una forma de calcularlo es usar estos valores pasados como variables dependientes en un modelo predictivo. Este enfoque es usado en un modelo basado en el suavizado exponencial: \[\begin{align} \hat{x}_{n+1} &= \alpha x_n+\alpha(1-\alpha)x_{n-1}+\alpha(1-\alpha)^2x_{n-2}+\cdots\\ &= \alpha x_n + (1-\alpha)\hat{x}_n \end{align}\]

donde $-1<\alpha<1$, asegurando que el factor $\alpha(1-\alpha)^i$ vaya disminuyendo a medida de que $i$ crece. Nótese que estas ponderaciones forman una serie geométrica. Para cualquier $\alpha$ dado, esta puede ser usada para predecir un valor $\hat{x}_t$ correspondiente a un valor $x_t$ observado, para $t=1,2,\cdots,n$. El error residual en la predicción está dado por: \[\begin{equation} \epsilon_t=x_t-\hat{x}_t \end{equation}\]

Y un buen estimador para $\alpha$ se obtiene minimizando la suma cuadrada de los errores: \[\begin{equation} \label{eq_SSE} SSE = \sum_{t=1}^n \epsilon_t^2 = \epsilon_1^2 + \epsilon_2^2+\cdots+\epsilon_n^2 \end{equation}\]

El método de Holt-Winters generaliza la ecuación anterior y simplifica la suavización al considerar sólo dos exponenetes suavizantes, permitiendo un término de tendencia cambiante utilizando el parámetro $\beta$ y un término de estacionariedad cambiante utilizando el parámetro $\gamma$. Para una serie $\{x_t\}$ con periodo $p$, tendencia $m_t$, y componente estacional $s_t$, el procedimiento Holt-Winters actualiza la tendencia y el efecto estacional así pronosticando el valor de $x_{n+1}$: \[\begin{equation} \hat{x}_{n+1} = m_n+b_n+s_{n+1-p} \end{equation}\] donde $b_n$ es el cambio estimado de la tendencia a tiempo $n$, entonces $m_n+b_n$ es la tendencia esperada a tiempo $n+1$, y $s_{n+1-p}$ es el efecto estacional estimado usando el valor anterior.

\[m_n =\alpha(x_n-s_{n-p})+(1-\alpha)(m_{n-1}+b_{n-1})\] \[b_n=\beta(m_n-m_{n-1})+(1-\beta)b_{n-1}\] \[s_n =\gamma(x_n-m_{n})+(1-\gamma)s_{n-p}\]

En , la función se usa para estimar los parámetros del modelo de Holt-Winters por minimización de suma de cuadrados del error (SSE ). El modelo de suavización exponencial también puede ajustarse con esta función, con los argumentos de y igualados a cero.

Entonces aplicamos el método Holt-Winters a los datos y hacemos una predicción

plot(CONELDO.ts)

Ahora se pasa la serie temporal a HoltWinter y se trazan los datos ajustados.

hw <- HoltWinters(CONELDO.ts)
plot(hw)

hw

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = CONELDO.ts)
## 
## Smoothing parameters:
##  alpha: 0.5744525
##  beta : 0.008084536
##  gamma: 0.3111858
## 
## Coefficients:
##           [,1]
## a   4094.67342
## b      7.55825
## s1  -308.73880
## s2  -437.69387
## s3  -324.13317
## s4   -45.48406
## s5   164.79579
## s6   374.71300
## s7   449.90018
## s8   486.25465
## s9   273.08280
## s10   16.82568
## s11 -284.28051
## s12 -229.77167

A continuación, calculamos el consumo de energia electrica en 24 meses con un intervalo de confianza de 0,95 y trazamos el pronóstico junto con los valores reales y ajustados.

forecast <- predict(hw, n.ahead = 24, prediction.interval = T, level = 0.95)
plot(hw, forecast)

plot(fitted(hw))

Además de esto R tiene incorporado el módelo de Holt winters multiplicativo. Por defecto R realiza el aditivo, el cual se puede escribir como:

hw_aditivo <- HoltWinters(CONELDO.ts, seasonal = "additive")
plot(hw)

Si se desea el módelo múltiplicativo se realiza:

hw_multiplicativo <- HoltWinters(CONELDO.ts, seasonal = "multiplicative")
plot(hw_multiplicativo)

Además en la minimización de las diferencias se puede dar valores por defecto a $\alpha, \beta, \gamma$ siempre que estos esten contenidos entre 0 y 1

parametros=HoltWinters(CONELDO.ts, beta=0, gamma=0)
plot(parametros)

parametros

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = CONELDO.ts, beta = 0, gamma = 0)
## 
## Smoothing parameters:
##  alpha: 0.4723123
##  beta : 0
##  gamma: 0
## 
## Coefficients:
##            [,1]
## a   4029.959013
## b      2.807255
## s1  -303.788194
## s2  -420.871528
## s3  -350.121528
## s4  -101.121528
## s5   176.211806
## s6   301.170139
## s7   473.336806
## s8   490.295139
## s9   211.211806
## s10   23.128472
## s11 -260.663194
## s12 -238.788194

Ruido Blanco

Como se mencionó anteriormente un residual $\varepsilon_t$ es la diferencia entre el valor real y el valor estimado en el t-ésimo momento, es decir: \[\varepsilon_t=x_t-\hat{x}_t\] En donde $x_t$ es el valor observado de la serie en el tiempo $t$ y $\hat{x}_t$ es el valor estimado en el mismo periodo.

De esta manera los residuales de la serie $\{\varepsilon_t: \ t=1, 2, \ldots, n\}$ es un ruido blanco si la variables $\varepsilon_1, \varepsilon_2,\ldots, \varepsilon_n$ son independientes e identicamente distribuidas con media 0 y varianza $\sigma^2$ y $cor(\varepsilon_i, \varepsilon_j)=0$ para toda $i\neq j$.

De igual forma, se el ruido blanco Gaussiano, la cual es el caso particular de un ruido blanco que sigue una distribución normal, $\varepsilon_t \sim N(0, \sigma^2) \forall t =1,\ldots, n$. De esta menera se concluye que para una serie de tiempo de tamaño $n$ tiene $n$ ruidos blancos, $\left(\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_n\right)$ de forma que cumple que:

\[E[ \varepsilon_t ]=\mu_\varepsilon=0 \]
\[Cov(\varepsilon_t, \varepsilon_{t+k})=\gamma_k=\begin{cases} \sigma^2 \quad si \quad k=0 \\ 0 \quad si \quad k\neq0 \end{cases} \]

Por lo que esto ocaciona que

\[\rho_k=\frac{Cov(\varepsilon_t, \varepsilon_{t+k})}{Cov(\varepsilon_t, \varepsilon_{t})}=\begin{cases} 1 \quad si \quad k=0 \\ 0 \quad si \quad k\neq0 \end{cases} \]

Estos resultados son teóricos, en la práctica puede resultar que la autocorrelación no sea exactamente 0, cuando $k\neq0$. Sin embargo, se esperaría que menos del 5% de las observaciones $k\neq0$ fuera significativamente igual a 0 con un nivel de confianza del 95%, ($\alpha=0.05$).

Para ejemplificar observe que:

set.seed(20201)
ruido_blanco=rnorm(100, 0, 5)
plot(ruido_blanco, type='l', main='Comportamiento')
abline(h=0, col='red')

hist(ruido_blanco)

acf(ruido_blanco)

Caminatas aleatorias

Sea $x_t$ una serie de tiempo. Entonces $x_t$ es una caminata aleatoria si: \[x_t=x_{t-1}+\varepsilon_t\]

Sustituyendo $x_{t-1}$ por el valor de la caminata se se puede concluir en una forma recursirva de la forma \[x_t=\varepsilon_t+\varepsilon_{t-1}+\varepsilon_{t-2}+\ldots\] Dado que la serie es finita y para $x_{1}=x_{0}+\varepsilon_1$ se tiene que por convención $x_{0}=0$ se tiene que:

\[x_t=\varepsilon_t+\varepsilon_{t-1}+\varepsilon_{t-2}+\ldots+\varepsilon_2+\varepsilon_1\]

A esta sustitución se le conoce como sustitución hacia atrás y es usado para modelos que cumplen con las propiedades de ruido blanco. Para ejemplificar veamos:

##COn el ruido blanco, se genera una caminata aleatoria usando el sustitución hacía atrás
##X(t)=X(t-1)+e(t)
##X(t)=X(t-2)+e(t-1)+e(t)
##Así hasta 
##X(t)=e(1)+e(2)+...+e(t-1)+e(t)
##
set.seed(20201)
RuidoBlanco=rnorm(355, mean=0, sd=25)
CaminataAleatoria=RuidoBlanco[1]
  for (i in 2:355) {
    CaminataAleatoria[i] = sum(RuidoBlanco[1:i])
  }
  
plot(CaminataAleatoria, type = "l")

Vemos que esta grafica es equivalente a la obtenida por la fórmula: \[x_t=x_{t-1}+\varepsilon_t\]

x=e=RuidoBlanco
for (t in 2:355) x[t] = x[t-1] + e[t]
plot(x, type="l")

Estacionariedad de la caminata aleatoria

Debido a los propiedades del ruido blanco la caminata aleaotoria tiene las siguientes propiedades: \[E[x_t]=\mu_x=0\] \[Cov[x_t, x_{t+k}]=t\sigma^2\] Por lo que la autocorrelación resultante es \[\rho_k(t)=\frac{cov(x_t, x_{t+k})}{\sqrt{Var(x_t)Var(x_{t+k})}}\] \[\rho_k(t)=\frac{t\sigma^2}{\sqrt{t\sigma^2(t+k)\sigma^2}}\]

\[\rho_k(t)=\frac{1}{\sqrt{1+k/t}}\]

Se observa que de manera general, para cualquier caminata aleatoria se tiene que $\rho_0=1$. Para el siguiente valor $\rho_k<\rho_{k+1}$, es decir la autocorrelación siempre es decreciente ya que depende exclusivamente de cada cardinalidad de la serie de tiempo

acf(CaminataAleatoria)

Operador Diferencia hacia atrás

El opereador de diferencia hacia atrás, denotado como $\phi$, en el curso haremos de el como operador lag, el cual es definido como: \[\phi x_t= x_{t-1}\] Se puede generalizar al operador lag como:

\[\phi^n x_{t}= x_{t-n}\]

Para entender mejor su uso se puede reescribir la caminata aleatoria como:

\[\begin{align} xt&=\phi x_t + \varepsilon\\ \varepsilon&=xt-\phi x_t \\ \therefore xt&= \frac{\varepsilon}{1-\phi} \end{align}\]

De igual se puede reescribir la caminata aleatoria con sistitución hacia atrás, debe entenderse que $\phi$ no esta definida sólo para la serie de tiempo, por lo que puede asumirse que :

\[\phi \varepsilon_t= \varepsilon_{t-1}\] \[\phi^2 \varepsilon_t= \varepsilon_{t-2}\] \[\vdots\] \[\phi^n \varepsilon_t= \varepsilon_{t-n}\]

Se observa así que :

\[x_t=\varepsilon_{t}+\varepsilon_{t-1}+\ldots+\varepsilon_1\] \[x_t=\varepsilon_{t}+\phi \varepsilon_{t}+\ldots+\phi^{t-1}\varepsilon_1\] \[x_t=(1+\phi+\phi^2+\ldots+\phi^{t-1})\varepsilon_t\]

Estacionariedad

Como habíamos mencionado anteriormente, una propiedad importante que se desea en una serie de tiempo es que sea estacionaria, ya que esto nos induce que las observaciones entre sí son linealmente independientes en el tiempo, en el caso de una caminata aleatoria por su construcción se tiene las observaciones estan muy relacionadas la una de la otra, ya su simulación es recursiva.

Un método para convertir una serie no estacionaria en una estacionaria es calculando las diferencias entre un periodo determinado. Sea el operador de diferencias $\nabla$, en este curso lo llamaremos operador dorito, de tal forma que:

\[\nabla x_t = x_t - x_{t-1}\]

Cabe destacar que al realizar esta operación se pierde un valor de la serie de tiempo, por lo que es necesario aplicar el operador dorito sólo cuando se tengan un tamaño considerable de datos, el valor perdido sería cuando $t=1$ ya que para \[\nabla x_1 = x_1 - x_{0},\] $x_0$ no esta definido.

Una caminata aleatoria es un proceso estocastico $x_t$ sí la diferencia de esta caminata aleatoria con su valor anterior es ruido blanco, es decir: \[\nabla x_t= \varepsilon_t\]

Realicemos esto en R.

##Calculo de diferencias con la función diff
proceso_estocastico=diff(CaminataAleatoria)
plot(proceso_estocastico, type='l')
abline(h=0, col='red')

acf(proceso_estocastico)

El operador $\nabla$ puede ser definida de una manera más generalizada como: \[\nabla^n x_t = x_t-x_{t-n}\]

Esta generalización es muy ocupada para hacer estacionaria a una serie de tiempo, ya que es usual que $n$ sea ajustada por la frecuencia, aunque no es una regla, es el valor con el debes intentar primerameramente. Por ejemplo, la serie de consumo de energia electrica en México, tiene tedencia y periodicidad de 12 datos por año, por lo que se puede realizar una estandarización.

##Calculo de diferencias
CONELDO.estaci=diff(CONELDO.ts, lag=12)


par(mfrow=c(2,2))
plot.ts(CONELDO.ts, main='Serie original')
acf(CONELDO.ts, main='AutoCorr. original')
plot.ts(CONELDO.estaci, main='Serie Transformada')
acf(CONELDO.estaci, main='AutoCorr. Transformada')

par(mfrow=c(1,1))

Modelos autorregresivos

La serie $x_t$ es un proceso autorregresivo de orden $p$, abreviado como $AR(p)$ si: \[x_t=\alpha_1 x_{t-1}+\alpha_2 x_{t-2}+\ldots+\alpha_p x_{t-p}+\varepsilon_t\] Donde $\varepsilon$ es ruido blanco, y $\alpha$ son paramétros tal que $\alpha_p \neq 0$. Usando la proipedad del Lag $\phi$ se puede expresar $AR(P)$ como: \[(1-\alpha_1 \phi - \alpha_2 \phi^2-\ldots- \alpha_p \phi^p)x_t=\varepsilon_t\]

El módelo autoregresivo es importante ya que explica varias soluciones a trevés de esta generalización. Por ejemplo:

La caminata aleatoria es un caso especial de un modelo $AR(1 )$ con $\alpha_1=1$ y $\alpha_i=0$ para toda $i=1,\ldots, p$ \[x_t=x_{t-1}+\varepsilon_t\]
La suavización exponencial es un caso especial cuando $\alpha_i=\alpha(1-\alpha)^i$ para toda $i=1, \ldots,$ \[x_t=\alpha(1-\alpha) x_{t-1}+\alpha(1-\alpha)^2 x_{t-2}+\ldots + \varepsilon_t\]
Un módelo de regresión puede ser usado usando $x_{t_1} como variable explicativa de $x_{t}$ sin ordendada al origen, si $\alpha_1=\beta$ y $\alpha_i=0$ para toda $i=1,\ldots, p$ \[x_t=\beta x_{t-1}+\varepsilon_t\] Para predicir una observación se denotará como $\hat{x}_t=$ el cual corresponde a \[\hat{x}_t=\alpha_1 x_{t-1}+\alpha_2 x_{t-2}+\ldots+\alpha_p x_{t-p}+\varepsilon_t\]

Los paramátros son estimados por mínimos cuadrados ordinaros, buscando mínimizar la suma de las distancias entre una valor real y una predicción.

\[SE=\sum_{i=1}^n (x_i-\hat{x}_i)^2\]

Para ejemplificar vemos el caso particular de AR(1)

Modelo AR(1)

EL modelo $AR(1)$ autorregresivo con un parámetro corresponde a una particularidad del modelo $AR(P)$, el cual se define como:

\[AR(1): x_t=\alpha x_{t-1}+\varepsilon_t\]

Este modelo puede ser expresado como: \[x_t-\alpha x_{t-1}=\varepsilon_t\] \[x_t(1-\alpha \phi)=\varepsilon_t\]

\[x_t=(1-\alpha \phi)^{-1}\varepsilon_t\] Otra manera de expresar el modelo AR(1) sería recursivamente como: \[x_t=\alpha x_{t-1}+\varepsilon_t\] \[x_t=\alpha (\alpha x_{t-2}+\varepsilon_{t-1})+\varepsilon_t\] Así sucesivamente hasta

\[x_t=\sum_{i=0}^{t-1} \alpha^i \varepsilon_{t-i}\]

De igual forma, se puede examinar que el proceso autorregresivo es una serie estacionaria ya que su esperanza y varianza son constantes en el tiempo:

\[E[x_t]=0\]
\[Cov(x_t, x_{t+k}))=\gamma_k=\alpha^k \sigma^2 \] es un preceso estacionario, ya que su esperanza y varianza son constantes en el tiempo:

El modelo AR(1) es muy parecido a un módelo de regresión lineal simple sin intercepto, el método para estimar $\alpha$ es por mínimos cuadrados ordinarios, en el que se busca mínimizar las distancias de los valores observados y los reales.

Si se minimiza la suma de cuadrados de la diferencia de los valores observados y los estimados

\[\left(\sum_{i=1}^{n} e_i^2\right),\]
Entonces se tiene como estimador de $\alpha$ a:

\[\hat{ \alpha }=\frac{\sum_{ i=1 }^{n}x_ix_{i-1}}{\sum_{i=1}^{n}x_i^2}\]

Demostración:

Se busca minimizar $\sum_{i=1}^{n} e_i^2$ por ello:

\[S(\alpha)=\sum_{i=1}^{n} e_i^2.\]

\[S(\alpha)=\sum_{i=1}^{n} (x_{i}-\hat{\alpha} x_{i-1})^2\] Derivando respecto a $\hat{\alpha}$ \[\frac{\partial S(\alpha)}{\partial \hat{\alpha}}=-2\sum_{i=1}^{n} (x_{i}-\hat{\alpha} x_{i-1})x_i\] Igualando la derivada a 0 para hallar el punto crítico.

\[-2\sum_{i=1}^{n} (x_{i}-\hat{\alpha} x_{i-1})x_i=0\] \[\sum_{i=1}^{n} x_ix_{i}-\hat{\alpha} x_{i-1}^2=0\] \[ \therefore \hat{\alpha}=\frac{\sum_{ i=1 }^{n}x_ix_{i}}{\sum_{i=1}^{n}x_i^2}\]

Volviendo a derivar para obtener si es máximo o mínimo.

\[\frac{\partial^2S(\alpha)}{\partial \hat{\alpha}^2}=(-2\sum_{i=1}^{n}(x_{i}x_i-\hat{\alpha} x_{i-1}^2))'\] \[\therefore \frac{\partial^2S(\alpha)}{\partial \hat{\alpha}^2} = 2\sum_{i=1}^{n}(x_{i}-\hat{\alpha} x_{i-1})^{-1}x_i^2>0.\]

Por lo tanto es un mínimo, por consiguiente, el estimador de $\alpha$ que minimiza la suma de cuadrados de los residuales es

\[ \hat{\alpha} =\frac{\sum_{ i=1 }^{n}x_ix_{i-1}}{\sum_{i=1}^{n}x_i^2}\] Las propiedades anteriormente descrita nos lleva a que la autocorrelación esta dado por:

\[\gamma_k=\alpha^k \forall k\ge 0\]

Con $\alpha \in (0,1)$. Lo que implica que la autocorrelación decae más rápido a 0 para valores cercanos a 0. Tiene sentido, ya que valores menores implica una relación no tan directa con la observación inmediata anterior, y valores más altos implica darle más peso a la observación anterior, por lo que la seríe estaría más autocorrelacionada, para ejemplificar, simulemos varios modelos autocorrelacionados

set.seed(3)
x = e = rnorm(100)
for (t in 2:100) x[t] = 0.01*x[t-1] + e[t]
par(mfrow=c(3,2))
plot(x, type='l', main='Simulación alpha=0.01')
acf(x, main='alpha=0.01')


for (t in 2:100) x[t] = 0.5*x[t-1] + e[t]
plot(x, type='l', main='Simulación alpha=0.5')
acf(x, main='alpha=0.5')


for (t in 2:100) x[t] = .99*x[t-1] + e[t]
plot(x, type='l', main='Simulación alpha=0.99')
acf(x, main='alpha=0.99')

par(mfrow=c(1,1))

Debido a esta problematica se creo otra medida conocida como la autocorrelación parcial (pacf). La autocorrelación parcial es el resultado de eliminar el efecto de cualquier correlación debido a términos anteriores. Por ejemplo, la autocorrelación parcial de un proceso AR (1) será cero para todos los rezagos mayores que 1. En general, la autocorrelación parcial en el rezago k es el coeficiente k de un modelo AR (k) ajustado; si el proceso subyacente es AR (p), entonces los coeficientes k serán cero para todos los k> p. Por lo tanto, un proceso AR (p) tiene un correlograma de autocorrelaciones parciales que “corta” a cero después del retraso p. Por lo tanto, una gráfica de las autocorrelaciones parciales puede ser útil al determinar el orden de Un proceso AR subyacente

set.seed(3)
x = e = rnorm(100)
for (t in 2:100) x[t] = 0.01*x[t-1] + e[t]
par(mfrow=c(3,2))
plot(x, type='l', main='Simulación alpha=0.01')
pacf(x, main='alpha=0.01')


for (t in 2:100) x[t] = 0.5*x[t-1] + e[t]
plot(x, type='l', main='Simulación alpha=0.5')
pacf(x, main='alpha=0.5')


for (t in 2:100) x[t] = .99*x[t-1] + e[t]
plot(x, type='l', main='Simulación alpha=0.99')
pacf(x, main='alpha=0.99')

par(mfrow=c(1,1))

La autocorrelación parcial será muy útil para estimar el grado del modelo autorregresivo, ya que las líneas que salen, indican las observaciones que estan altemente autocorrelacionadas entre si, por lo que al elegir ese grado se dará mejores proyecciones. Para ejemplificar veamos:

set.seed(3)
x = e = rnorm(100)
for (t in 2:100) x[t] = 0.5*x[t-1] + e[t]
par(mfrow=c(2,1))
pacf(x, main='alpha=0.01')


for (t in 3:100) x[t] = 0.2*x[t-1]+ 0.5*x[t-2] + e[t]
pacf(x, main='alpha=0.5')

par(mfrow=c(1,1))

Una vez idetentificado el orden entonces, se puede estimar el parametro alpha. La forma automatizada es simplemente poner el metódo a usar en este caso, se usará el método de mínimos cuadrados ordinarios

set.seed(3)
x = e = rnorm(100)
for (t in 2:100) x[t] = 0.5*x[t-1] + e[t]


x.ar = ar(x, method="mle")
x.ar$order; x.ar$ar

## [1] 1

## [1] 0.5203917

for (t in 3:100) x[t] = 0.2*x[t-1]+ 0.5*x[t-2] + e[t]
x.ar = ar(x, method="mle")

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

x.ar$order; x.ar$ar

## [1] 2

## [1] 0.2356326 0.4854777

Este proceso es automático, sin embargo, si se desea un valor en especifico para AR(p) se puede establecer como:

set.seed(3)
x = e = rnorm(100)
for (t in 3:100) x[t] = 0.2*x[t-1]+ 0.5*x[t-2] + e[t]
x=ts(x, start = 1)
x.ar=ar (x,  order.max = 1)
x.ar$order; x.ar$ar

## [1] 1

## [1] 0.4541088

Para ver los valores predichos puede verse que:

set.seed(3)
x = e = rnorm(100)
for (t in 3:100) x[t] = 0.8*x[t-1]+ .2*x[t-2] + e[t]
x=ts(x, start = 1)
x.ar=ar (x)
x.ar$order; x.ar$ar

## [1] 1

## [1] 0.9176344

prediccion=predict(x.ar, x, n.ahead = 25)
prediccion$pred

## Time Series:
## Start = 101 
## End = 125 
## Frequency = 1 
##  [1]  1.26286114  1.04596762  0.84693865  0.66430281  0.49670988  0.34292084
##  [7]  0.20179871  0.07230019 -0.04653211 -0.15557672 -0.25563981 -0.34746115
## [13] -0.43171957 -0.50903800 -0.57998805 -0.64509427 -0.70483797 -0.75966085
## [19] -0.80996821 -0.85613198 -0.89849344 -0.93736578 -0.97303637 -1.00576894
## [25] -1.03580547

plot(x, col='blue')
lines(prediccion$pred, col='red')

plot(x, col='blue', type='l', xlim=c(80,110))
lines(prediccion$pred, col='red')

Aplicación con series reales

La serie del consumo de energia electrica en México vemos que no estacionaria, por lo que aplicamos diferencia hacia atrás, buscando obtener una serie estacionaria

plot(CONELDO.ts)

CONELDO.dorito=diff(CONELDO.ts, lag=12)
plot(CONELDO.dorito)

Una vez estacionario el módelo se hace el ajuste de series de tiempo, y se observa la autocorrelación parcial para saber que parámetros elegir

pacf(CONELDO.dorito)

Como se ve que salen dos líneas entonces probamos intentar realizar un modelo AR(2).

CONELDO.ar=ar(CONELDO.dorito, order.max = 2)
CONELDO.ar

## 
## Call:
## ar(x = CONELDO.dorito, order.max = 2)
## 
## Coefficients:
##      1  
## 0.4193  
## 
## Order selected 1  sigma^2 estimated as  11681

Se proyecta la información todo un año

prediccion=predict(CONELDO.ar, CONELDO.dorito, n.ahead = 12)

plot(CONELDO.dorito, col='blue', type='l', xlim=c(2010, 2019 ))
lines(prediccion$pred, col='red')

library(forecast)

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

pronostico_coneldo=forecast(CONELDO.ar, 11, level=95)
plot(pronostico_coneldo, main="pronostico CONELDO")

matriz.pronostico_coneldo=data.frame(pronostico_coneldo$mean, pronostico_coneldo$lower, pronostico_coneldo$uppe)
matriz.pronostico_coneldo

##    pronostico_coneldo.mean       X95.   X95..1
## 1                319.02177  107.18754 530.8560
## 2                189.05031  -40.65125 418.7519
## 3                134.55430  -98.14670 367.2553
## 4                111.70454 -121.51978 344.9289
## 5                102.12382 -131.19239 335.4400
## 6                 98.10669 -135.22566 331.4391
## 7                 96.42235 -136.91285 329.7575
## 8                 95.71611 -137.61959 329.0518
## 9                 95.41999 -137.91579 328.7558
## 10                95.29583 -138.03997 328.6316
## 11                95.24377 -138.09203 328.5796

escritura_coneldo=pronostico_coneldo$mean
diferencias=c(CONELDO.dorito,escritura_coneldo)
serie_orginal=diffinv(diferencias, lag=12, differences = 1,CONELDO.ts[1:12] )

CONELDO.ts_add=ts(serie_orginal,start =2009,frequency = 12)
CONELDO.ts_add

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2009 2861.000 2833.000 2686.000 2798.000 3123.000 3232.000 3377.000 3545.000
## 2010 2825.000 2769.000 2655.000 2735.000 2996.000 3282.000 3594.000 3544.000
## 2011 2998.000 3062.000 2907.000 3012.000 3349.000 3537.000 3688.000 3632.000
## 2012 2999.000 2990.000 2983.000 3068.000 3279.000 3440.000 3636.000 3670.000
## 2013 3126.000 3120.000 2905.000 2992.000 3272.000 3524.000 3732.000 3807.000
## 2014 3230.000 3176.000 3005.000 3119.000 3382.000 3543.000 3820.000 3985.000
## 2015 3229.000 3212.000 3043.000 3188.000 3483.000 3678.000 3896.000 4142.000
## 2016 3415.000 3330.000 3063.000 3284.000 3677.000 3964.000 4211.000 4296.000
## 2017 3371.000 3331.000 3297.000 3396.000 3753.000 3967.000 4322.000 4341.000
## 2018 4000.000 3650.022 3486.050 3530.554 3864.705 4069.124 4420.107 4437.422
##           Sep      Oct      Nov      Dec
## 2009 3558.000 3275.000 3079.000 2792.000
## 2010 3633.000 3422.000 3220.000 2859.000
## 2011 3792.000 3608.000 3274.000 3049.000
## 2012 3835.000 3605.000 3436.000 3126.000
## 2013 3793.000 3631.000 3481.000 3203.000
## 2014 3950.000 3743.000 3512.000 3211.000
## 2015 4190.000 4063.000 3847.000 3482.000
## 2016 4310.000 4142.000 3858.000 3591.000
## 2017 4398.000 4213.000 3867.000 3509.000
## 2018 4493.716 4308.420 3962.296 3604.244

require(xts)
plot(as.xts(CONELDO.ts_add), major.ticks = 'years', main='Consumo de energia domestica')

Módelo Medias Móviles

Un módelo de medias móviles de orden $q$, denotado como $MA(q)$ estad definido como: \[x_t=\varepsilon_t+\theta_1 \varepsilon_{t-1}+\ldots+\theta_q \varepsilon_{t-q}, \] donde $\varepsilon_j\sim N(0, \sigma^2), \ j>0$

El módelo puede escribirse con los operadores anteriormente vistos como:

\[x_t=(1+\theta_1 \psi+\theta_2\psi^2+\ldots+\theta_q\psi^q)\varepsilon, \]

La suma de uno hasta $\theta_q\psi^q$ puede escribirse como:

\[x_t=\theta(\psi)\varepsilon, \] donde $\theta(\psi)$ es el operador de media móvil. Al igual que como sucedio en el módelo AR, se particulariza escencialmente en el módelo AM(1)

Media Móvil de grado 1

El módelo de media móvil de grado 1, es el caso más sencillo, es denotado como $MA(1)$, lo que denota: \[x_t=\varepsilon_t+\theta \varepsilon_{t-1}\]

Este módelo cumple con ser estacionario, ya que: \[E[x_t]=0\] \[Var[x_t]=(1+\theta^2)\sigma^2\] \[Cov[x_t, x_{t-1}]=\theta \sigma^2\] \[Cov[x_i, x_{i+k}]=0 , \forall k>1\]

Resumiendo la función de autocovarianzas es:

\[\gamma_k=\begin{cases} (1+\theta^2)\sigma^2 & si \ k=0 \\ \theta \sigma^2 & si \ k=1 \\ 0 & si \ k>0\end{cases}\]

En este caso debido a que las autocovarianzas no son dependientes en el tiempo se establece que el método para estimar el grado de la media móvil, son aquellas observaciones autocorrelacionadas entre si menos 1, debido a que la primera autocorrelación siempre es 1.

Simulación

Para ejemplificar un modelo AM(1), supongo que desea estimar un modelo de medias móviles, con $\theta=0.4$ y $\varepsilon_i\sim N(0,1)$, de esta forma, la serie de tiempo es la siguiente:

\[x_t=\varepsilon_t+0.4 \varepsilon_{t-1}.\] La simulación de 100 valores es el siguiente:

Ejemplo

library(ggplot2)
set.seed(3)
x = e = rnorm(100)
for (t in 2:100) x[t] = 0.4*e[t-1] + e[t]
plot(x, type="l", main="Simulación de series de tiempo")
abline(h=0, col="red")

serie_tiempo=data.frame(time=1:100,serie=x)
ggplot(serie_tiempo, aes(x=time, y=serie))+
  geom_line()+geom_line( color="steelblue") + 
  geom_hline(yintercept=0, linetype="dashed")+
  theme(axis.text.x=element_text(angle=60, hjust=1))+
  labs(x="Tiempo", y="Serie de tiempo", title = "Simulación serie de tiempo")

Se verifican las propiedades, como que:

\[E[x_t]=0\] \[\gamma_k=\begin{cases} (1+0.4^2)\sigma^2=1.16 & si \ k=0 \\ \theta \sigma^2=0.4 & si \ k=1 \\ 0 & si \ k>0\end{cases}\]

esperanza=mean(x)
esperanza

## [1] 0.0162869

varianza=sd(x)
varianza

## [1] 0.9250104

rho1=cov(x[1:99], x[2:100])
rho1

## [1] 0.3020268

rho2=cov(x[1:98], x[3:100])
rho2

## [1] 0.01673267

rho50=cov(x[1:50], x[51:100])
rho50

## [1] -0.01715121

La función de autocorrelacion se define como \[\rho_k =\frac{\gamma(k)}{\gamma(0)}=\begin{cases} 1 & si \ k=0 \\ \frac{\theta}{1+\theta^2}& si \ k=1 \\ 0 & si \ k>1\end{cases} \]

Se obseva que $0\le \rho_1 \le 0.5$, así realizando las autocorrelaciones.

autocorrelacion=acf(x)

autocorrelacion

## 
## Autocorrelations of series 'x', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.349  0.019  0.070  0.078  0.089 -0.052 -0.107  0.035  0.042 -0.003 
##     11     12     13     14     15     16     17     18     19     20 
## -0.029 -0.043 -0.004 -0.072 -0.109 -0.096 -0.057 -0.159 -0.204 -0.088

Para estimar, el módelo se usará el método de máxima verosímilitud, el cual se aproximará a través de R.

Modelo=arima(x, order=c(0,0,1), include.mean = FALSE)
Modelo

## 
## Call:
## arima(x = x, order = c(0, 0, 1), include.mean = FALSE)
## 
## Coefficients:
##          ma1
##       0.4192
## s.e.  0.0954
## 
## sigma^2 estimated as 0.7245:  log likelihood = -125.88,  aic = 255.76

El módelo estimado por R, es el siguiente: \[x_t=\varepsilon_t+0.4192\varepsilon_{t-1}\]

Lo cual es un valor muy acercado al valor propuesto el que fue propuesto. Finalmente se proyectará la información 15 periodos hacia adelante.

library(forecast)
pronostico=forecast(Modelo, 15, level=95)
plot(pronostico, main="Pronóstico")

Modelo AM(2)

El modelo AM(2) se puede escribir como: \[x_t=\varepsilon_t+\theta_1 \varepsilon_{t-1}+\theta_2\varepsilon_{t-2}\]

La simulación de 100 valores es el siguiente:

Ejemplo

library(ggplot2)
set.seed(8)
x = e = rnorm(100)
for (t in 3:100) x[t] = 0.2*e[t-1] + .6*e[t-2]+ e[t]
plot(x, type="l", main="Simulación de series de tiempo")
abline(h=0, col="red")

serie_tiempo=data.frame(time=1:100,serie=x)
ggplot(serie_tiempo, aes(x=time, y=serie))+
  geom_line()+geom_line( color="steelblue") + 
  geom_hline(yintercept=0, linetype="dashed")+
  theme(axis.text.x=element_text(angle=60, hjust=1))+
  labs(x="Tiempo", y="Serie de tiempo", title = "Simulación serie de tiempo")

La forma en que esta construido el módelo se tiene que:

\[E[x_t]=0\] \[\gamma_k=\begin{cases} (1+\theta_1^2+\theta_2^2)\sigma^2 & si \ k=0 \\ \theta_1 \sigma^2(1+\theta_2) & si \ k=1 \\ \theta_2\sigma^2 & si \ k=2\\ 0 & si \ k>2\end{cases}\] Si $\theta_1=0.2$, $\theta=0.4$, se tiene que: \[\gamma_k=\begin{cases} (1+0.2^2+0.4^2)=1.2 & si \ k=0 \\ 0.2 (1+0.4)=.28 & si \ k=1 \\ 0.4 & si \ k=2\\ 0 & si \ k>2\end{cases}\]

esperanza=mean(x)
esperanza

## [1] -0.1634149

varianza=sd(x)
varianza

## [1] 1.101226

rho1=cov(x[1:99], x[2:100])
rho1

## [1] -0.02147419

rho2=cov(x[1:98], x[3:100])
rho2

## [1] 0.3148371

rho3=cov(x[1:97], x[4:100])
rho3

## [1] -0.05640771

rho50=cov(x[1:50], x[51:100])
rho50

## [1] 0.1423626

La función de autocorrelacion se define como \[\rho_k =\frac{\gamma(k)}{\gamma(0)}= \begin{cases} 1 & si \ k=0 \\ \frac{\theta_1(1+\theta_2)}{1+\theta^2_1+\theta^2_2} & si \ k=1 \\ \frac{\theta_2}{1+\theta^2_1+\theta^2_2}& si \ k=2 \\ 0 & si \ k>2 \end{cases}\]

Se obseva que $0\le \rho_1 \le 2/3$, y $0\le \rho_2 \le 1/3$ así realizando las autocorrelaciones.

autocorrelacion=acf(x)

autocorrelacion

## 
## Autocorrelations of series 'x', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.018  0.254 -0.045 -0.077  0.163  0.090 -0.005 -0.006 -0.349 -0.078 
##     11     12     13     14     15     16     17     18     19     20 
## -0.126  0.069  0.050 -0.041 -0.144 -0.016  0.014  0.159  0.071  0.011

Para estimar, el módelo se usará el método de máxima verosímilitud, el cual se aproximará a través de R.

Modelo=arima(x, order=c(0,0,2), include.mean = FALSE)
Modelo

## 
## Call:
## arima(x = x, order = c(0, 0, 2), include.mean = FALSE)
## 
## Coefficients:
##          ma1     ma2
##       0.0659  0.4017
## s.e.  0.0897  0.1041
## 
## sigma^2 estimated as 1.094:  log likelihood = -146.57,  aic = 299.13

El módelo estimado por R, es el siguiente: \[x_t=\varepsilon_t+0.0659\varepsilon_{t-1}+0.4017\varepsilon_{t-2}\]

Lo cual es un valor muy acercado al valor propuesto el que fue propuesto. Finalmente se proyectará la información 15 periodos hacia adelante.

library(forecast)
pronostico=forecast(Modelo, 15, level=95)
plot(pronostico, main="Pronóstico")

Ejercicio serie real

Para ejemplificar continuemos con el consumo de energía eléctrica en México (CONELDO). Puntualmente trabajamos con el módelo estacionario que vimos anteriormente (diferencia anual)

plot(CONELDO.dorito)

acf(CONELDO.dorito)

Vemos que de la gráfica de 4 líneas sobresalen, sin embargo, debido a que el primero siempre es 1, se sugiere que el modelo a estimar es una MA(3), así:

Modelo=arima(CONELDO.dorito, order = c(0,0,3))
Modelo

## 
## Call:
## arima(x = CONELDO.dorito, order = c(0, 0, 3))
## 
## Coefficients:
##          ma1     ma2     ma3  intercept
##       0.5634  0.1807  0.1152    99.3324
## s.e.  0.1434  0.1234  0.1494    19.4557
## 
## sigma^2 estimated as 10722:  log likelihood = -587.89,  aic = 1185.78

El modelo obtenido es el siguiente \[x_t=99.3324+0.5634\varepsilon_{t-1}+0.1807 \varepsilon_{t-2}+0.1152 \varepsilon_{t-3}+\varepsilon_t\]

Al igual que en el módelo pasado buscamos estimar todo el año 2018. Así

library(forecast)
pronostico_coneldo=forecast(Modelo, 11, level=95)
plot(pronostico_coneldo, main="pronostico CONELDO")

matriz.pronostico_coneldo=data.frame(pronostico_coneldo$mean, pronostico_coneldo$lower, pronostico_coneldo$uppe)
matriz.pronostico_coneldo

##    pronostico_coneldo.mean       X95.   X95..1
## 1                418.78112  215.83405 621.7282
## 2                196.08183  -36.85969 429.0234
## 3                171.38256  -64.42803 407.1932
## 4                 99.33239 -137.63421 336.2990
## 5                 99.33239 -137.63421 336.2990
## 6                 99.33239 -137.63421 336.2990
## 7                 99.33239 -137.63421 336.2990
## 8                 99.33239 -137.63421 336.2990
## 9                 99.33239 -137.63421 336.2990
## 10                99.33239 -137.63421 336.2990
## 11                99.33239 -137.63421 336.2990

escritura_coneldo=pronostico_coneldo$mean
diferencias=c(CONELDO.dorito,escritura_coneldo)
serie_orginal=diffinv(diferencias, lag=12, differences = 1,CONELDO.ts[1:12] )

CONELDO.ts_add=ts(serie_orginal,start =2009,frequency = 12)
CONELDO.ts_add

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2009 2861.000 2833.000 2686.000 2798.000 3123.000 3232.000 3377.000 3545.000
## 2010 2825.000 2769.000 2655.000 2735.000 2996.000 3282.000 3594.000 3544.000
## 2011 2998.000 3062.000 2907.000 3012.000 3349.000 3537.000 3688.000 3632.000
## 2012 2999.000 2990.000 2983.000 3068.000 3279.000 3440.000 3636.000 3670.000
## 2013 3126.000 3120.000 2905.000 2992.000 3272.000 3524.000 3732.000 3807.000
## 2014 3230.000 3176.000 3005.000 3119.000 3382.000 3543.000 3820.000 3985.000
## 2015 3229.000 3212.000 3043.000 3188.000 3483.000 3678.000 3896.000 4142.000
## 2016 3415.000 3330.000 3063.000 3284.000 3677.000 3964.000 4211.000 4296.000
## 2017 3371.000 3331.000 3297.000 3396.000 3753.000 3967.000 4322.000 4341.000
## 2018 4000.000 3749.781 3493.082 3567.383 3852.332 4066.332 4421.332 4440.332
##           Sep      Oct      Nov      Dec
## 2009 3558.000 3275.000 3079.000 2792.000
## 2010 3633.000 3422.000 3220.000 2859.000
## 2011 3792.000 3608.000 3274.000 3049.000
## 2012 3835.000 3605.000 3436.000 3126.000
## 2013 3793.000 3631.000 3481.000 3203.000
## 2014 3950.000 3743.000 3512.000 3211.000
## 2015 4190.000 4063.000 3847.000 3482.000
## 2016 4310.000 4142.000 3858.000 3591.000
## 2017 4398.000 4213.000 3867.000 3509.000
## 2018 4497.332 4312.332 3966.332 3608.332

require(xts)
plot(as.xts(CONELDO.ts_add), major.ticks = 'years', main='Consumo de energia domestica')

Arima

Un poceso Arima es una generalización de un módelo AM y uno AR, ya que se involucra cada uno de esos procesos, en su construcción, de ahí es donde se deriva el nombre, ARIMA (“AR” de un proceso autoregresivo, “MA” de un proceso de medias móviles e “i” indica el grado $\nabla$ que deben calcularse). De esta manera se define a la serie de tiempo de tamaño $n$ como un proceso ARIMA(p, d, q) como: \[x_t=-(\nabla^d x_t - x_t)+\alpha_1 \nabla^d x_{t-1}+\alpha_2 \nabla^d x_{t-2}+\ldots+\alpha_p \nabla^d x_{t-p}+\varepsilon_t+\theta_1\varepsilon_1+\ldots+\theta_q\varepsilon_q\] El procedimiento aunque pareciera complicado no lo es, ya que la anterior formula se puede resumir si primero se calcula el grado del operador dorito a tal grado que la serie sea estacionaria. Es decir primero se calcula $\nabla^d x_t=x_t-x_{t-d}$, una vez encontrado aquel grado que haga estacionario a la serie de tiempo, ie, se tenga el claro patrón de un ruido blanco, se puede expresar la serie como: \[x_t=\alpha_1 x_{t-1}+\alpha_2 x_{t-2}+\ldots+\alpha_p x_{t-p}+\varepsilon_t+\theta_1\varepsilon_1+\ldots+\theta_q\varepsilon_q\]

Lo que es equivalente a \[x_t=AR(p)+MA(q)-\varepsilon_t\]

Series de tiempo

Omar Rodríguez Torres

11/03/2021

Introducción

Notación

Ejemplo 2

Ejemplo pasajeros de avión

Correlación y Correlograma

Covarianza

Estacionariedad y doble estacionaridad

Descomposición de modelos

Modelo aditivo

Estamación de tendencia y estacionalidad

Pronósticos

Ruido Blanco

Caminatas aleatorias

Estacionariedad de la caminata aleatoria

Operador Diferencia hacia atrás

Estacionariedad

Modelos autorregresivos

Modelo AR(1)

Aplicación con series reales

Módelo Medias Móviles

Media Móvil de grado 1

Simulación

Modelo AM(2)

Ejercicio serie real

Arima