Estimación de la tendencia y la componente estacional: método S1

Recordemos que en este caso estimamos la tendencia calculando el promedio por año:

\[ \hat{\mu_{j}}=\frac{1}{12}\sum_{k=1}^{12}x_{j,k}, \ \ \ \ \ j=1,\ldots,r \]

es decir, fijamos un año y el promedio corre sobre el ciclo estacional \(k\). De forma similar estimamos la componente estacional como

\[ \hat{s_{k}}=\frac{1}{r}\sum_{j=1}^{12}(x_{j,k}-\hat{\mu_{j}}), \ \ \ \ \ k=1,\ldots,r \]

es decir, fijamos el ciclo estacional \(k\) y corremos el promedio sobre los años, donde en este promedio quitamos la parte de la tendencia para cada respectivo año. Finalmente, para la parte aleatoria quitaremos la parte d ela tendencia y de la componente estacional.

Ejemplo

Utilizaremos los datos de USAccDeaths (muertes accidentales en EU)

# Cargamos la librería necesaria
library(forecast)

# guardamos los datos
xt = USAccDeaths

# vemos la gráfica de la ST
plot(xt)

notamos que el ciclo se repite cada año y los datos abarcan 6 años (\(j=1,\ldots,6\)). Así, utilizando el método \(S1\) estimamos la componente de la tendencia (\(\hat{\mu_{j}}\)) como sigue:

  • Tomando \(j=1\) realizaremos la suma de los datos de la ST sobre los meses y dividimos dicha suma entre 12
sum( xt[1:12] ) / 12
  • Para \(j=2\) haremos lo mismo, pero ahora nos interesan los datos del segundo año, o lo que es lo mismo, nos interesan los datos de la segunda fila:
sum( xt[12 + 1:12] ) / 12

donde recordemos que la primer fila de la ST es de la forma \(x_{1}, x_{2},\ldots,x_{12}\) y la segunda fila \(x_{13}, x_{14},\ldots,x_{24}\), de ahí que hayamos colocado xt[12 + 1:12] para acceder a estos datos de la ST. De forma general

# creamos un vector vacío para que almacene los valores
mj = c()

# usamos un for para recorrer todos los años de los datos de la ST
for(j in 1:6){
  mj[j] = sum( xt[12*(j-1)+1:12] ) / 12
}
mj
## [1] 9651.750 8718.500 8588.583 8395.083 8576.833 8802.000

Hacemos lo mismo para estimar la componente estacional

# vector vacío
sk = c()

# usamos un for para recorrer todos los meses
for(k in 1:12){
  sk[k] = sum( xt[(0:5)*12+k]-mj ) / 6
}

donde xt[(0:5)*12+k]-mj nos permite acceder a los datos de la ST con índice \(1,13,25,37,49,61\) para \(k=1\) (Enero), etcétera.

Finalmente estimamos la parte aleatoria

aleatorio = xt - rep(mj, each = 12) - rep(sk, times = 6)

donde, como los datos de la ST son de 6 años con periodo de 12 meses, entonces podemos pensar a los datos como una tabla de 6 filas por 12 columnas. De tal manera, a la ST xt le restaremos entrada por entrada los valores obtenido es mj y sk.

Podemos ver graficamente la parte aleatoria obtenida y la gráfica de la ST

par(mfrow=c(2,1))
plot(xt)
plot(aleatorio)

asimismo, podemos graficar, en azul, la estimación de la componente de la tendencia junto con la estimación de la componente estacional

# compenente estacional y de tendencia (estimación)
compon = rep(mj, each = 12) + rep(sk, times = 6) 

# lo adecuamos a los datos de la ST
compon = ts(compon, start=start(xt), frequency=12)

# graficamos
plot(xt)
lines(compon,col="blue")

Estimación de la tendencia y la componente estacional: método S2

Recordemos que para este caso debemos identificar primero el periodo de la parte estacional, de modo que en este caso \(d=12\); además, dado que la ST abarca 6 años, tenemos que \(k=6\). De tal manera, \(n=k\cdot d=6\cdot 12=72\). En efecto

length(xt)
## [1] 72

Ahora, dado que \(d\) es par, tenemos entonces $12 =2q $, de donde \(q=6\). Luego

d = 12
k = 6
n = 72
q = 6

tenemos entonces que el estimador de la tendencia está dado por:

# Definimos una ST sólo con entradas "NA" de acuerdo al número de 
# entradas de la ST original
mt = ts(rep(NA, n), start = start(xt), frequency = frequency(xt))

# Realizamos la suma
for(t in (q+1):(n-q)){
  mt[t] = (0.5 * xt[t - q] + sum( xt[(t-q+1):(t+q-1)] ) + 0.5 * xt[t + q]) / d  
}

donde recordemos que \(q+1\leq t\leq n-q\) y además la fórmula para estimar la tendencia es

\[ \hat{m_{t}}=\frac{0.5 x_{t-q}+x_{t-q+1}+\cdots+x_{t+q-1}+0.5x_{t+q}}{d} \]

dado que en las primeras \(q\) entradas de la ST no estamos obteniendo esta estimación de la tendnecia, y dado que tampoco estamos realizando esta estimación para los últimos datos (i.e \(t>n-q\)), hacemos que en ambos casos la tendencia sea constante

mt[1:q] = mt[q+1]
mt[(n-q+1):n] = mt[n-q]

Lo que sigue es quitar la tendencia a la serie original (\(z_{t}=x_{t}-\hat{\mu_{t}}\))

zt = xt - mt

Después, a la ST auxiliar \(z_{t}\) crearemos un ciclo promedio (\(w_{k}\), \(k=1,\ldots, 12\)) que estimará la parte estacional:

wk  = rep(NA, 12)

for(k in 1:12){
  wk[k] = sum(zt[(0:5)*12 + k], na.rm = TRUE) / 6
}

donde para \(k=1\) tendríamos que promediar los datos \(z_{1}, z_{13}, z_{25}, z_{37}, z_{49}, z_{61}\), etcétera.

Continuamos ajustando el ciclo obtenido en la serie como

\[ \hat{S_{k}}=w_{k}-\frac{\sum_{i=1}^{d}w_{i}}{d}, \ \ \ k\in\{1,\ldots,12\} \]

# Vec_ciclo:

sk  = rep(NA, 12)

for(k in 1:12){
  sk[k] = wk[k] - sum(wk) / d
}

Finalmente repetimos el vector \(Vec_{ciclo}\) tantos ciclos tengamos (a lo cual denominaremos \(S_{t}\)), en este caso repetiremos el vector sk 6 veces:

# Adaptamos los datos a una serie de tiempo
St = ts(rep(sk, times=6), start=start(xt), frequency=frequency(xt))

Podemos calcular la parte aleatoria como \(y_{t}=z_{t}-S_{t}\)

aleatorio = zt-St

Finalmente vemos los resultados:

  1. La ST original junto con la estimación de la tendencia y la estimación del ciclo
# Componente de la tendencia más componente estacional
compon = mt + St

# Lo adaptamos a una ST
compon = ts(compon, start=start(xt), frequency = 12)

# Graficamos
plot(xt)
lines(compon,col="red")

y finalmente podemos observar la tendencia estimada, el ciclo estimado y la parte aleatoria (respectivamente):

par(mfrow=c(3,1))
plot(mt, col="blue")
plot(St, col="blue")
plot(aleatorio)