Descomposición clásica y suavizamiento exponencial

Pregunta 1.

Grafique los datos. Describa lo que observe de su información (varianza contante o no constante, tendencia, ciclos estacionales, periodicidad de los ciclos).

Solución: Comenzamos por cargar y ver la serie de tiempo:

# Cargamos las librerias que ocuparemos
library(forecast)
library(astsa)
library(nortest)
library(lmtest)

# Cargamos los datos
data <- scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")

# Adecuamos los datos a una serie de tiempo
data_TS <- ts(data, frequency = 12, start = c(1946, 1))

# Graficamos la serie de tiempo utilizando plot()
plot(data_TS)

Podemos observar a primera instancia que la varianza es constante, pues a simple vista no vemos un crecimiento o decrecimiento entorno a la dispersión; además vemos que la tendencia es no lineal.

Notamos además la presencia de un ciclo estacional con periodicidad de cada 12 meses, lo cual podemos sustentar valiéndonos de las correlaciones y correlaciones parciales

par(mfrow=c(1,2))
acf(data_TS)
pacf(data_TS)

Pregunta 2.

Use distintos métodos de descomposición de las series para obtener sus componentes (tendencia y ciclo estacional), en específico use los siguientes

\(a)\) Ajuste de curvas (Modelo de regresión). Realice un pronóstico de 3 años futuros.

\(b)\) Promedios móviles, filtros lineales o suavizamientos exponenciales. Realice un pronóstico de 3 años futuros

Inciso a)

Solución: \(a)\) Por simple inspección de la gráfica de la serie de tiempo podemos intuir que la tendencia sigue un cierto comportamiento polinomial. Por ello

t = time(data_TS) - 1945 
M <- factor(cycle(data_TS))

t2 = t ** 2
t3 = t ** 3
t4 = t ** 4
t5 = t ** 5

# Realizamos el ajuste del modelo
reg <- lm( data_TS ~ 0 + t + t2 + t3 + t4 + t5 + M , na.action = NULL )

# Vemos alguna informacion sobre este ajuste
summary(reg)
## 
## Call:
## lm(formula = data_TS ~ 0 + t + t2 + t3 + t4 + t5 + M, na.action = NULL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01321 -0.42624 -0.00468  0.40833  1.84254 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## t   -8.431e+00  1.002e+00  -8.411 2.85e-14 ***
## t2   2.448e+00  3.423e-01   7.153 3.41e-11 ***
## t3  -3.071e-01  5.127e-02  -5.989 1.48e-08 ***
## t4   1.819e-02  3.472e-03   5.238 5.37e-07 ***
## t5  -4.136e-04  8.693e-05  -4.757 4.54e-06 ***
## M1   3.151e+01  9.925e-01  31.750  < 2e-16 ***
## M2   3.002e+01  9.956e-01  30.153  < 2e-16 ***
## M3   3.303e+01  9.985e-01  33.079  < 2e-16 ***
## M4   3.136e+01  1.001e+00  31.324  < 2e-16 ***
## M5   3.245e+01  1.004e+00  32.325  < 2e-16 ***
## M6   3.200e+01  1.006e+00  31.801  < 2e-16 ***
## M7   3.359e+01  1.008e+00  33.313  < 2e-16 ***
## M8   3.340e+01  1.011e+00  33.050  < 2e-16 ***
## M9   3.290e+01  1.013e+00  32.494  < 2e-16 ***
## M10  3.296e+01  1.015e+00  32.484  < 2e-16 ***
## M11  3.106e+01  1.017e+00  30.550  < 2e-16 ***
## M12  3.180e+01  1.019e+00  31.215  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7179 on 151 degrees of freedom
## Multiple R-squared:  0.9993, Adjusted R-squared:  0.9992 
## F-statistic: 1.213e+04 on 17 and 151 DF,  p-value: < 2.2e-16
anova(reg)
## Analysis of Variance Table
## 
## Response: data_TS
##            Df Sum Sq Mean Sq   F value    Pr(>F)    
## t           1  90106   90106 174821.33 < 2.2e-16 ***
## t2          1  10460   10460  20293.40 < 2.2e-16 ***
## t3          1   2969    2969   5761.03 < 2.2e-16 ***
## t4          1   1444    1444   2801.19 < 2.2e-16 ***
## t5          1    636     636   1234.89 < 2.2e-16 ***
## M          12    704      59    113.76 < 2.2e-16 ***
## Residuals 151     78       1                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

donde se puede decir que el modelo es bastante bueno pues \(R^{2}=0.9993\) y además el \(ECM\) de los residuales es igual a 1. Si agregamos más variables a la regresión, el \(ECM\) aumenta y la \(R^{2}\) disminuye, por lo que se ha decido trabajar con el modelo presentando anteriormente. Luego, graficamos el ajuste del modelo de la serie de tiempo junto con la gráfica de los datos observados de la misma:

plot(data_TS, type = "o")
lines(fitted(reg), col = 2, lwd = 3)
legend(x = "topleft", c("Original", "ajuste"), fill = c("black", "red"))

De lo anterior podemos obtener la tendencia y el ciclo estacional de la serie de tiempo.

Tendencia

Realizaremos otro ajuste del modelo, pero esta vez excluiremos a la variable categórica M referente al ciclo estacional de la serie de tiempo, así

# Realizamos el ajuste del modelo sin la variable M
reg_trend <- lm( data_TS ~ 0 + t + t2 + t3 + t4 + t5 , na.action = NULL )

# graficamos la serie de tiempo
plot(data_TS, type = "o")

# graficamos la linea correspondiente a la tendencia de la serie de tiempo
lines(fitted(reg_trend), col = 2, lwd = 3)
legend(x = "topleft", c("Serie de tiempo", "tendencia"), fill = c("black", "red"))

Ciclo estacional

Para obtener el ciclo estacional de la serie de tiempo aplicaremos primero la función diff() a la serie para quitar la tendencia y así visualizar mejor el ciclo estacional. Después, realizaremos otro ajuste del modelo quitando las variables t:

# Realizamos el ajuste del modelo sin las variables t
reg_cicl <- lm( data_TS ~ 0 + M , na.action = NULL )

# Graficamos la serie de tiempo
plot(data_TS, type = "o")

# graficamos las lineas correspondientes al ciclo de la serie de tiempo
lines(fitted(reg_cicl), col = 2, lwd = 3)
legend(x = "topleft", c("Serie de tiempo", "ciclo estacional"), fill = c("black", "red"))

Finalmente, para la parte de predicción escribimos

# pronostico
pred <- predict(fitted(reg), h=36)

# graficamos
par(mfrow=c(1,2))
plot(pred)

Inciso b)

\(b)\) Para este inciso utilizaremos el método \(S2\), por ello identificamos

length(data_TS)
## [1] 168
d = 12
k = 14
n = 168
q = 6

de donde \(12=2\cdot q\), por lo que \(q=6\). Tenemos entonces que el estimador de la tendencia está dado por:

# Definimos una ST solo con entradas "NA" de acuerdo al numero de 
# entradas de la ST original
mt = ts(rep(NA, n), start = start(data_TS), frequency = frequency(data_TS))

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

dado que en las primeras \(q\) entradas de la ST no estamos obteniendo esta estimación de la tendencia, y dado que tampoco estamos realizando esta estimación para los últimos datos (i.e \(t>n-q\)), hacemos que en ambos casos que 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 = data_TS - 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(k2 in 1:12){
  wk[k2] = sum(zt[(0:5)*12 + k2], na.rm = TRUE) / d
}

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(k2 in 1:12){
  sk[k2] = wk[k2] - 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 14 veces:

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

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 mas componente estacional
compon = mt + St

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

# Graficamos
plot(data_TS)
lines(compon,col="red")
legend(x = "topleft", c("Original", "ajuste"), fill = c("black", "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)

Pregunta 3.

Aplique transformaciones y diferencias para hacer varianza constante y quitar tendencias y ciclos estacionales.

Solución: Transformamos usando la función \(log\), y quitamos los ciclos estacionales y la tendencia usando la función Diff().

# Transformamos
data_TS_log <- log(data_TS)

# Vemos la grafica de los datos transformados
plot(data_TS_log)

con lo cual vemos más notoriamente que la varianza se comporta de forma constante, además, parecería que el ciclo continúa. Procedemos entonces a aplicar la función Diff()

plot(diff(data_TS_log))

con lo cual hemos conseguido quitar la tendencia. Nuevamente aplicamos la función Diff() para quitar el ciclo estacional

tsdisplay((diff(diff(data_TS_log, 12))))

donde, con base en ACF y PACF, vemos que ya no hay un ciclo estacional.

Vemos las varianzas para ver sí nuestra transformación y diferenciación es útil. Para ello, observamos primero las varianzas con los datos de la serie de tiempo original

# original
var(data_TS)
## [1] 5.376791
# quitando tendencia
var(diff(data_TS))
## [1] 2.279232
# quitando ciclo estacional
var(diff(diff(data_TS,12)))
## [1] 0.929462

y después las varianzas con la serie de tiempo transformada

# datos trasformados
var(data_TS_log)
## [1] 0.008707438
# quitando tendencia
var(diff(data_TS_log))
## [1] 0.003747315
# quitando ciclo estacional
var(diff(diff(data_TS_log,12)))
## [1] 0.001500268

Pregunta 4.

Use el método de Holt-Winter para el ajuste de la curva y predicción de los datos de 3 años futuros.

Solución: Comenzamos aplicando un suavizamiento exponencial a la serie del tiempo, además de visualizar al gráfica de la serie de tiempo junto con el ajuste de la curva proveniente del método Holt-Winters

# Suavizamiento Exponencial
data_TS_exp <- HoltWinters(data_TS)

# Vemos la grafica
plot(data_TS_exp)
legend(x="topleft", legend=c("Original", "Ajuste: Holt-Winters"), fill=c("black", "red"))

Para la parte de predicción utilizaremos el paquete forecast y la función forecast.HoltWinters():

# Cargamos la libreria
library("forecast")

# Definimos una variable para la prediccion
data_TS_exp_pred <- forecast:::forecast.HoltWinters(data_TS_exp, h=36)

# Graficamos
plot(data_TS_exp_pred)
legend(x="topleft", legend=c("Original", "Predicción a 3 años"), fill=c("black", "blue"))