Gráficos de Series de Tiempo

Una forma de ver la evolución de datos de series de tiempo es generando un gráfico de series de tiempo, en este se puede ver la tendencia que ha ido siguiendo cierta variable a través del tiempo. La función plot() permite hacer gráficos para variables de series de tiempo en R.

plot(gnp.ts, ylab="Producto Nacional Bruto", xlab="Trimestres")

Componentes de una Serie de Tiempo

Una serie de tiempo es un dato que se registra a través del tiempo. Este tipo de datos tiene a su vez componentes propios, suponiendo que son 4 componentes que actúan al mismo tiempo, son los que contribuyen a formar las características de dicha serie, estos cuatro componentes son la tendencia, la estacionalidad, la ciclicidad y la variación irregular o aleatoriedad.

Descomposición de Series de Tiempo

Antes de hablar de descomposición de una serie de tiempo se debe entender que existen modelos que explican el comportamiento de una serie de tiempo.

gnp.decomp = decompose(gnp.ts, type = "mult")

tend = gnp.decomp$trend
estac = gnp.decomp$seasonal
aleat = gnp.decomp$random

plot(gnp.decomp)

Procesos Estocásticos Estacionarios

Sea la serie \(W_t\). La estacionariedad sirve para limitar o restringir lavheterogeneidad a través del tiempo de un proceso estocástico.

1. De forma estricta la estacionariedad hace referencia a que las funciones de distribución conjunta de subconjuntos de un número cualquiera 𝑛 de una variable no varía por traslaciones en el tiempo para cualquier horizonte ℎ:

\[ f(w_{t1}, w_{t2},...,w_{tn})=f(w_{t1+h},w_{t2+h},...,w_{tn+h}) \forall n,h \]

2. En un sentido más amplio un proceso es estacionario si cumple las siguientes condiciones:

• \(E(w_t) = \mu, \forall t\)
• \(0 < Var(W_t) = \sigma^2 < \infty, \forall t\)
• \(Cov(W_t, W_{t+k}) = \gamma (k),k = \pm1, \pm2,..., \forall t\)

Cuando el proceso de análisis tiene una distribución normal, la estacionariedad en sentido estricto coincidirá con la de sentido amplio. La estacionariedad se obtiene diferenciando la serie y aplicando logaritmos (Trasformación Box-Cox). A partir de las correlaciones se va a poder tomar decisiones sobre cuál es el tipo de modelo estadístico más adecuado para dicha serie de tiempo. Cuando se impone la estacionariedad, se reduce número de parámetros, pero todavía queda un número elevado de parámetros por \(\gamma_k\), se necesitará una condición más, que es la ergodicidad (consistencia). Cuando se quiere hacer pruebas de hipótesis y análisis de los parámetros, se va a dar con que a medida que aumenta el número de observaciones, aumentará el número de parámetros desconocidos. La condición de ergodicidad será la siguiente: La correlación serial entre variables de una serie disminuye a medida que nos alejamos en el tiempo. \[ \lim_{x\rightarrow\infty} \gamma_k = 0 \]

Proceso de Ruido Blanco

Un proceso de ruido blanco o White noise, es el modelo más simple de series temporales. Este modelo trata de una serie puramente aleatoria y se representa por 𝑎𝑡 .Una serie será de ruido blanco cuando:

• La esperanza es constante, , es igual a cero.
• \(Cov(a_t, a_{t+k})= 0 \text{ para todo } k \neq 0\)

Este proceso se trata básicamente de uno en el cual todas sus variables son independientes. Un ejemplo de un proceso de ruido blanco son los números de una lotería:

• Cada número es independiente del anterior. • No hay dependencia entre el pasado y el futuro.

Es decir, la variable tiempo no influye, ya que no hay dependencia, esto es totalmente impredecible. Un modelo de series temporales adopta la siguiente expresión: \[ y_t=f(pasado) + a_t = f(y_{t-1},t_{t-2},...) + a_t \] En 𝑡 − 1 la función es conocida, pero 𝑎𝑡 será desconocida, ya que es una variable aleatoria que no se puede predecir.

set.seed(1)
w = rnorm(100)

plot(w, type="l")

Proceso Random Walk

El proceso random walk viene dado por la siguiente expresión: \[ y_t = y_{t-1} + a_t \] Donde \(a_t\) es un proceso de ruido blanco, es decir, entra en el proceso como shock aleatorio que se incorpora en cada momento. Este proceso tiene raíz unitaria, (coeficiente de \(y_{t-1}\)), por lo tanto, la serie muestra un perfil o una tendencia evolutiva, por lo que se entenderá que no es estacionaria. Como ejemplo de este tipo de procesos, podemos notar a diferentes variables económicas en mercado eficientes, en este tipo de mercados se entiende que con un elevado número de agentes con información completa los valores presentes de las series se adaptarán a la información disponible, es decir el precio de cierto bien en el periodo \(t-1\) será el mismo que el de futuro.

y <- w <- rnorm(1000)

for (t in 2:1000) y[t] <- y[t - 1] + w[t]

plot(y, type = "l")

La Autocorrelación

acf(gnp.ts)$acf[2]

## [1] 0.9759739

acf(gnp.ts, type = c("covariance"))$acf[2]

## [1] 2799951

El Correlograma

Por defecto, a función acf(), produce gráfico de \(r_k\) contra \(k\) que es llamada correlograma. Por ejemplo, si escribimos la sintaxis acf(gnp.ts)obtendremos el correlograma para dicha serie.

acf(gnp.ts)

De este gráfico se puede entender que:

El eje \(x\) da el rezago \(k\) y el eje \(y\) da la autocorrelación \(r_k\) en cada rezago. Si \(p_k = 0\), la distribución muestral de \(r_k\) es aproximadamente normalm con una media de \(-1/n\) y una varianza de \(1/n\). Las líneas del correlogama se grafican en: \[ -\frac{1}{n}\pm\frac{2}{\sqrt{n}} \] Si \(r_k\) cae afuera de esas líneas, se tiene evidencia contra la hipótesis nula \(p_k = 0\)al 5% de significancia. Sin embargo se debe tener cuidado al interpretar la hipótesis nula. En primer lugar, si \(p_k\) es igual a 0 en el rezago \(k\) , se espera que el 5% de los estimados caiga fuera de las líneas. En segundo lugar, os \(r_k\) están correlacionados, por lo tanto, si uno cae fuera de las líneas, los coeficientes cercanos están más probables a ser estadísticamente singificativos.

Correlograma de un ruido blanco

El correlograma para un proceso de ruido blanco tiene la siguiente figura:

Correlograma

Correlograma para un proceso de random walk

El correlograma para un proceso de paseo aleatorio o random walk tiene la siguiente forma:

Correlograma de un Ruido Blanco

Correlogama de un Proceso AR(1)

La función de autocorrelación de un proceso AR(1) viene dada por:

\[ p_k=\alpha^k (k\geq0) \] Donde \(|\alpha| < 1\). Así, el correlograma decrece a 0 o más rápidamente para los \(\alpha\) pequeños. El siguiente ejemplo muestra dos correlogramas para valores positivos y negativos de \(\alpha\), respectivamente. La sintaxis siguiente mostrará el siguiente gráfico:

rho <- function(k, alpha) alpha^k
layout(1:2)
plot(0:10, rho(0:10, 0.7),type = "b")
plot(0:10, rho(0:10, -0.7), type = "b")

Se puede ver que hay decaimiento rápido hacia 0.

Autocorrelación Parcial

La autocorrelación parcial en el rezago \(k\) es la correlación que resulta de quitar el efecto de cualquier correlación debido a los términos de los rezagos más cortos. Por ejemplo, la autocorrelación parcial de un proceso *AR(1) será 0 para todos los rezagos más grande que 1. En general, la autocorrelación parcial en el rezago \(k\) es el \(k\)vo coeficiente de un modelo AR(k), si el proceso subyacente es un AR(p)**, entonces los coeficientes \(\alpha_k\) serán cero para todo \(k > p\).

set.seed(1)
y <- w <- rnorm(100)
for (t in 2:100) y[t] <- 0.7 * y[t-1] + w[t]
layout(1:3)
plot(y, type = "l"); acf(y); pacf(y)

y.ar = ar(y, method = "ols")

y.ar$order

## [1] 1

y.ar$ar

## , , 1
## 
##           [,1]
## [1,] 0.6015796

Modelos Lineales

Un modelo para una serie de tiempo \((x_t:t = 1,...,n)\) es lineal si se puede expresar como: \[ x_t = \alpha_0+\alpha_1\mu_{1,t}+\alpha_2\mu_{2,t}+...+\alpha_m\mu_{m,t}+Z_t \] Donde 𝑢𝑖,𝑡 es el valor de la i-ésima variable predictora en el tiempo 𝑡, 𝑧𝑡 es el error en el tiempo 𝑡, y 𝛼0, 𝛼1, … , 𝛼𝑚 son parámetros del modelo, que pueden ser estimados por mínimos cuadrados. Un ejemplo del modelo lineal es la función polinomial de 𝑡 de orden 𝑝: \[ x_t = \alpha_0+\alpha_1t+\alpha_2t^2+...+\alpha_pt^p+Z_t \] Las variables predictoras pueden ser escritas como 𝑢𝑖,𝑡 = 𝑡𝑖 (𝑖 = 1, … , 𝑝). El término lineal es una referencia de la sumatoria de los parámetros del modelo, cada una por la variable predictora.

set.seed(1)
z = w = rnorm(100, sd=20)
for (t in 2:100) z[t] = 0.8 *z[t-1] + w[t]
Time = 1:100
x = 50+3*Time+z
plot(x, xlab = "time", type = "l")

Ajuste del Modelo

Los modelos lineales se suelen ajustar usualmente minimizando la suma de los cuadrados del error \(\sum Z_t^2 = \sum(x_t-\alpha_0-\alpha_1\mu_1-...-\alpha_m\mu_{m,t})^2\), que puede ser estimado en R con la función lm.

x.lm = lm(x~Time)
coef(x.lm)

## (Intercept)        Time 
##   58.551218    3.063275

sqrt(diag(vcov(x.lm)))

## (Intercept)        Time 
##  4.88006278  0.08389621

En la sintaxis indicada, se extraen los parámetros del modelo estimado con la función coef(). Se debe notar, como se esperaba, que los valores estimados son cercanos a los parámetros especificados, 50 para el intercepto y 3 para el coeficiente. Los errores estándar se extraen sacando la raíz cuadrada a los elementos de la diagonal de la matriz de varianzas y covarianzas (vcov). Después de ajustar el modelo de regresión, se deben considerar diversos gráficos de diagnóstico. En el caso de la regresión de series de tiempo, un diagnostico importante es el gráfico de correlograma de los residuos, haciendo uso de la función de autocorrelación y la función de autocorrelación parcial, de la siguiente manera:

acf(resid(x.lm))

pacf(resid(x.lm))

Como se esperaba, la serie de tiempo del residuo está autocorrelacionada, la figura de la autocorrelación parcial muestra que sólo en el primer rezago la autocorrelación parcial es significativa, lo que sugiere que la serie del residuo sigue un proceso 𝐴𝑅(1).

Modelos Lineales con Variables Estacionales

Ya que las series de tiempo son observaciones medidas de forma secuencial en el tiempo, suelen presentarse efectos estacionales en los datos.

Variables aditivas de indicadores estacionales

Suponga que una serie de tiempo contiene \(s\) estaciones. Por ejemplo, para una serie de tiempo medida para cada mes del calendario, \(s=12\), así como las series semestrales, con \(s=2\). Un modelo con indicador estacional para una serie de tiempo \((x_t:t=1,...,n)\) que contiene \(s\) estaciones y una tendencia \(m_t\) está dada por: \[ x_t = m_t + s_t + z_t \] Donde \(s_t=\beta_i\) cuando \(t\) cae en la i-ésima estación (t=1,…,n ; i = 1,…,s) y \(z_t\) es el error residual de la serie, que puede estar autocorrelacionada . Este modelo toma la misma forma que un modelo de descomposición aditivo. La ecuación anterior no muestra en \(m_t\) un término constante, es decir \(m_t\) puede ser un polinomio de orden \(p\) con parámetros \(\alpha_1,...,\alpha_p\). Esta ecuación es entonces equivalente a una tendencia polinómica en donde al término constante depende de la estación, por lo tanto, los parámetros \(\beta_1,...,\beta_s\) corresponden a los \(s\) posibles términos constantes. La ecuación puede ser escriba como: \[ x_t = m_t + \beta_{1+(t-1)mod \text{ s }} + Z_t \] Por ejemplo, con una serie de tiempo observada para cada mes del calendario comenzando con \(t=1\) como enero, un indicador estacional con tendencia de línea recta está dado por:

Ecuación

Los parámetros del modelo en esta ecuación pueden ser estimados por MCO o Mínimos Cuadrados Generalizados, tratando al término estacional \(S_t\) como un factor. En R, la función factor() puede ser aplicada para los índices estacionales extraídos usando la función cycle(). Con la base de datos de sesiones pasadas sobre el producto nacional bruto trimestral, se puede estimar el modelo de regresión estacional de la siguiente manera:

gnp.ts = ts(gnp96$gnp96, start=c(1967, 1), end=c(2002,2), frequency = 4)
gnpt = window(gnp.ts, start = 1990)

Se está indicando la extracción de una parte de la base total, que será desde 1990, luego se deben extraer los índices estacionales con la función cycle() y el tiempo con la función time(); luego de eso se hará la regresión del modelo indicando que los índices estacionales son una variable de factor, haciendo uso de la función del mismo nombre (factor).

Seas <- cycle(gnpt)
Time <- time(gnpt)
temp.lm <- lm(gnpt ~ 0 + Time + factor(Seas))

Los coeficientes se observarán usando la función coef(), los ceros usados en la fórmula indica que el modelo no tiene intercepto.

coef(temp.lm)

##          Time factor(Seas)1 factor(Seas)2 factor(Seas)3 factor(Seas)4 
##      268.6879  -528322.4276  -528333.8611  -528366.4147  -528348.1618

Modelos No Lineales

Los modelos lineales son aplicables en cierto rango. Sin embargo, para algunas series de tiempo puede ser más apropiado ajustar modelos no lineales de forma directa en lugar de usar logaritmos o usar aproximaciones polinómicas. Por ejemplo, si se conoce que una serie está derivada de un proceso no lineal, puede que esté basada en una ley determinística en la ciencia, es por eso que será mejor usar esta información en la formulación del modelo y ajustar un modelo no lineal directamente a los datos. Para estimar modelo no lineales en R, se debe hacer usar la función nls(), que hace referencia a la estimación por mínimos cuadrados no lineales. Se ha visto que se puede estabilizar la varianza aplicando logaritmos naturales a una serie. Sin embargo, usar logaritmos puede presentar dificultades cuando la serie presenta valores negativos, ya que el logaritmo de un valor negativo no está definido. Una manera de evitar dicho problema es añadiendo una constante a todos los términos de la serie, entonces, si \(X_t\) es una serie que contiene valores negativos, añadiendo \(c_o\) tal que \(c_o > max(x_t)\) y tomando los logaritmos se produce una serie transformada \((log(c_o+x_t))\) que está definida para todo \(t\). Un modelo lineal puede entonces ser ajustado para producir el modelo: \[ x_t = -c_o + e^{\alpha_0+\alpha_1t+z_t} \] Donde \(\alpha_0\) y \(\alpha_1\) son parámetros del modelo y \(Z_t\) es una serie residual que puede estar autocorrelacionada. La principal dificultad de la ecuación es que \(c_o\) debe ser estimado como otro parámetro del modelo, mientras que en la práctica el investigador suele elegir un valor arbitrario que satisface la condición \((c_o > max(-x_t))\). El modelo puede ser estimado por mínimos cuadrados no lineales, de la siguiente manera: \[ x_t = -c_o + e^{\alpha_0+\alpha_1t}+Z_t \] ## Simulación y ajuste de series no lineales

Como los modelos no lineales son generalmente ajustados cuando se conoce la función no lineal subyacente, se simulará una serie basada en la última ecuación, con \(c_o=0\) y se compararán los parámetros de la estimación usando la función nls().

set.seed(1)
w = rnorm(100, sd = 10)
z = rep(0 , 100)
for (t in 2:100) z[t] = 0.7 * z[t-1] + w[t]
Time <- 1:100
f <- function(x) exp(1 + 0.05 * x)
x <- f(Time) + z
plot(x, type = "l")
abline(0, 0)

El dibujo de la serie muestra un aparente crecimiento con tendencia exponencial que contiene valores negativos, por eso la transformación en logaritmos de la serie no puede ser usada y se necesita la especificación de un modelo no lineal. En R, un modelo no lineal es ajustado especificando una formula con los parámetros y los valores de entrada contenidos con la función list:

x.nls <- nls(x ~ exp(alp0 + alp1 * Time), start = list(alp0 = 0.1,
alp1 = 0.5))
summary(x.nls)$parameters

##        Estimate   Std. Error  t value     Pr(>|t|)
## alp0 1.17640120 0.0742954447 15.83410 9.204055e-29
## alp1 0.04828424 0.0008188462 58.96619 2.347543e-78

Los valores estimados de \(\alpha_0\) y \(\alpha_1\) están cerca a los valores subyacentes que se usaron para simular la data.

Pronósticos para Regresiones

Un pronóstico es una predicción para el futuro. En el contexto de regresión de series de tiempo, un pronóstico implica la extrapolación del modelo ajustado evaluando en él la función para nuevas series de tiempo. El problema principal con esta aproximación es que las tendencias que se presentan en las series ajustadas pueden cambiar en el futuro. Por lo tanto, es mejor pensar en un pronóstico de una regresión como un valor condicional esperado en tendencias pasadas que continúan en el futuro.

Seas <- cycle(gnpt)
Time <- time(gnpt)
temp.lm <- lm(gnpt ~ 0 + Time + factor(Seas))
#
new.t <- time(ts(start = 1990, end = c(2002, 2), fr = 4))
pron = ts(predict(temp.lm, new.t), st = 1990, fr=4)
ts.plot(gnp.ts, pron, lty = 1:2)

library(forecast)

## Warning: package 'forecast' was built under R version 4.0.2

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

mod <- tslm(gnpt ~ 0 + trend + season)
fc = forecast(mod, h = 8)

fc

##         Point Forecast    Lo 80     Hi 80    Lo 95    Hi 95
## 2002 Q3       9681.138 9407.565  9954.712 9257.500 10104.78
## 2002 Q4       9766.563 9492.990 10040.137 9342.924 10190.20
## 2003 Q1       9859.469 9585.423 10133.516 9435.098 10283.84
## 2003 Q2       9915.208 9641.161 10189.254 9490.837 10339.58
## 2003 Q3       9949.826 9673.699 10225.953 9522.233 10377.42
## 2003 Q4      10035.251 9759.124 10311.378 9607.658 10462.84
## 2004 Q1      10128.157 9851.380 10404.934 9699.558 10556.76
## 2004 Q2      10183.896 9907.119 10460.673 9755.296 10612.50

plot(fc, type = "l")

Ejemplo 1: Gráficos de Series de Tiempo

Se tiene la base de datos adjunta, que representa a la producción mensual de leche en libras por vaca, desde enero de 1962 hasta diciembre de 1975. En el siguiente ejemplo se verá cómo se declara cierta variable como serie de tiempo y cómo será el gráfico de la misma haciendo uso de R.

Luego de eso se deber declarar a R que la variable es de series de tiempo, para eso antes se modificarán los nombres de dichas variables.

View(milkpercow)
colnames(milkpercow) = c("mes", "milk")

Se declara a la variable milk como serie de tiempo:

milk.ts = ts(milkpercow$milk, start=c(1962, 1), end=c(1975,12), frequency = 12)

Una vez que se declara dicha variable indicando la frecuencia de la misma, se procede a hacer el gráfico de dicha serie de tiempo, para eso se usará la función plot().

plot(milk.ts, xlab = "Tiempo" , main = "Produccion mensual de leche")

Esta sintaxis le indicará a R que se quiere un gráfico de serie de tiempo donde el título será “Producción mensual de leche”, haciendo uso de la variable que se declaró como serie de tiempo.

Ejemplo 2: Descomposición de una serie de tiempo

Se tiene la base de datos adjunta, que representa al Producto Interno Bruto en millones de soles, desde el primer trimestre del año 2000 hasta el último del año 2015. En el siguiente ejemplo se realizará la descomposición de dicha serie para extraer sus componentes.

Luego de importar declaramos la variable como serie de tiempo:

pbi.ts = ts(pbi_trim$pib, start=c(2000, 1), end=c(2015,4) , frequency = 4)

Esto le habrá indicado a R que la variable pbi en la base de datos sigue una estructura de serie de tiempo que es trimestral. Una vez que se declara dicha variable indicando la frecuencia de la misma, se procede a hacer el gráfico de dicha serie de tiempo, para eso se usará la función plot(). Esta sintaxis le indicará a R que se quiere un gráfico de serie de tiempo donde el título será “Producto Interno Bruto (en millones de soles)”, haciendo uso de la variable que se declaró como serie de tiempo, también indicará que se quiere hacer un gráfico con líneas azules (opción col).

plot(pbi.ts, xlab = "Tiempo", main = "Producto Interno Bruto (en millones de soles)", col = "blue")

Una vez que se tiene el gráfico de la serie, se podrá observar cierto patrón de estacionalidad, para observar a más detalle los componentes de esta serie se deberá hacer la descomposición de la misma. Para eso se usará la función decompose(). Se asumirá que el modelo de esta serie es multiplicativo.

pbi.decomp = decompose(pbi.ts, type = "mult")

Esto creará un objeto que contendrá a cada uno de los componentes de la serie temporal. Se puede graficar cada uno de los componentes haciendo uso de la función plot().

plot(pbi.decomp, col = "brown")

Esta sintaxis indica a R que se quiere crear un gráfico del objeto que representa a la descomposición de la serie, se especifica el color de las líneas en el gráfico, que será marrón.

Ejemplo 4: Pronósticos de Series

El objetivo principal del ejemplo es mostrar el uso de la función predict(), que servirá para comparar el modelo estimado y los valores observados de la serie.

Primero declaramos a la serie cargada como una serie de tiempo:

library(tseries)

## Warning: package 'tseries' was built under R version 4.0.3

con.ts= ts(construct$pib, start=c(1950), end=c(2010))

Podemos hacer un gráfico de líneas de la serie en cuestión:

plot(con.ts)

Podemos observar que esta serie ha tenido un crecimiento exponencial, esto entonces puede ser explicado por un modelo no lineal, donde se deberá utilizar la estimación de mínimos cuadrados no lineales para encontrar los parámetros estimados. La forma del modelo a estimarse en R será la siguiente:

\[ const_t = e^{\alpha_0+\alpha_1t}+Z_t \] Se entiende que 𝑐0 = 0, y la sintaxis para estimar este tipo de modelos en R será la siguiente:

mod <- nls(con.ts ~ exp(alp0 + alp1 * time(con.ts)), start = list(alp0 = 0.1, alp1 = 0.5))

En la sintaxis se indica la forma del modelo y los valores de entrada, los parámetros podrán observarse haciendo de la función summary():

summary(mod)$parameters

##        Estimate   Std. Error   t value     Pr(>|t|)
## alp0 9.32364737 0.0406931731 229.12068 9.927974e-89
## alp1 0.02698092 0.0008582745  31.43623 1.498490e-38

Los resultados muestran que ambos parámetros estimados son estadísticamente significativos. Para estimar los pronósticos del modelo se debe crear una nueva estructura de datos, luego de esto creará una serie de tiempo con la variable pronósticos y finalmente se graficará los pronósticos.

new.t <- time(ts(start = 1950, end = c(2010)))
pron = ts(predict(mod, new.t), st = 1950)
ts.plot(con.ts, pron, lty = 1:2)