Cuando tenemos alguna explicación física plausible para una tendencia, podemos querer modelarla de manera determinística. Así, las tendencias determinísticas y la variación estacional regular pueden ser modeladas utilizando la regresión.
Comenzaremos a mirar los modelos lineales de tendencia y luego introducir modelos de regresión que den cuenta de la variación estacional utilizando variables indicadoras o armónicas. Asimismo, estudiaremos los modelos de regresión que incluyen variables explicatorias.
La regresión de series temporales difiere del análisis estándar de regresión dado que los residuales forman una serie temporal y por tanto tienden a estar correlacionados en la serie. Cuando esta correlación es positiva, los errores estándar estimados para los parámetros tenderán a ser menores que sus valores reales, conduciendo a una significación estadística erróneamente mayor. En estos casos, se utilizan mínimos cuadrados generalizados para obtener mejores estimaciones del error estándar y dar cuenta de la autocorrelación de las series residuales.
Nota: recordar que la transformación logarítmica estabiliza la varianza de la serie.
1. Modelos lineales:
Un modelo para la serie temporal {X(t), t=1,…,n} es lineal si puede expresarse como:
Nota: muchos modelos no-lineales pueden transformarse en modelos lineales. Por ejemplo, el modelo X(t)=exp(alfa_0+alfa_1*t+Z(t)) de la serie {X(t)} puede transformarse tomando logaritmos naturales para obtener un modelo lineal de la serie {Y(t)} mediante Y(t)=log(X(t)). Luego, para obtener predicciones sobre X(t), debemos utilizar la transformación inversa sobre Y(t): exp(Y(t)). (ver sección 5.10 para una discusión de este tema)
Nota: la diferenciación de la serie puede eliminar la tendencia estocástica y determinística: diff(X(t))=X(t)-X(t-1). Por ejemplo, si la tendencia subyacente es polinómica de orden m, necesitaremos una diferenciación de orden m para eliminar la tendencia.
Aplicación en R: Mínimos Cuadrados simples (OLS).
Aplicación en R: Mínimos Cuadrados Generalizados (GLS).
2. Modelo lineal con variables estacionales
En las series temporales pueden existir efectos estacionales debido especialmente a los ciclos anuales causados directa o indirectamente por el movimiento de la Tierra alrededor del Sol.
Ejemplo: para una serie temporal {X(t)} mensual (t=1, corresponde al mes de Enero), podemos estimar los parámetros del siguiente modelo mediante OLS o GLS, tratando al término estacional S(t) como factor para obtener los índices estacionales.
Es decir, el modelo es X(t)=alfa_1+S(t)+Z(t), donde S(t) es un factor con parámetros beta_1 (para t=1,13,…),…,beta_12 (para t=12,24,…).
Aplicación en R: Serie de Temperaturas (1970-2005)
Generalmente, se utilizan las funciones seno y coseno para construir la variación estacional del modelo estacional.
Ejemplo: Una ola seno con frecuencia f (ciclos por intervalos muestreado), de amplitud A y cambio de fase phi, se expresa como:
Aplicación en R: simularemos el modelo anterior de variables estacionales armónicas. Por ejemplo, para: X(t)= 0.1 + 0.005*t + 0.001*t^2 + sin(2*phi*t/12)+0.2*sin(4*phi*t/12) + 0.1*sin(8*phi*t/12) + 0.1*cos(8*phi*t/12) + W(t). con W(t) un ruido blanco gaussiano con desviación estándar .5.
Para series con valores negativos, podemos sumarle una constante a cada término de la serie, de tal manera que si {X(t)} es la serie que contiene algunos valores negativos, si le agregamos una constante c_0 tal que c_0>max{-X(t)} y tomamos logaritmos, obtenemos una serie transformada {log(c_0+X(t))} que está definida para todo t. Entonces, podemos ajustar un modelo lineal (e.g. tendencia lineal): x(t)=-c_0+exp(alfa_0+alfa_1*t+z(t)), donde alfa_0 y alfa-1 son los parámetros del modelo y {z(t)} es la serie residual (que puede estar autocorrelacionada). El problema es que el parámetro c_0 también debe ser estimado por el modelo, cuando la selección de c_0 suele hacerse tal que satisfaga la condición c_0>max(-X(t)). Si no existe evidencia para términos residuales multiplicativos, la constante c_0 será estimada con los demás parámetros del modelo, utilizando mínimos cuadrados no-lineales, ajustando el siguiente modelo: x(t)=-c_0+exp(alfa_0+alfa_1*t+z(t)).
3. Modelos no-lineales (e.g. modelo de Bass):
Para algunas series (e.g. cuando conocemos que el proceso subyacente es no-lineal) es más apropiado ajustar un modelo no-lineal directamente más que tomar logaritmos o utilizar aproximaciones polinómicas lineales.
Aplicación en R: simularemos una serie de datos no-lineal con c_0=0 y comparamos los parámetros estimados utilizando la función “nls” con aquellos parámetros conocidos de la función subyacente.
Comenzaremos a mirar los modelos lineales de tendencia y luego introducir modelos de regresión que den cuenta de la variación estacional utilizando variables indicadoras o armónicas. Asimismo, estudiaremos los modelos de regresión que incluyen variables explicatorias.
La regresión de series temporales difiere del análisis estándar de regresión dado que los residuales forman una serie temporal y por tanto tienden a estar correlacionados en la serie. Cuando esta correlación es positiva, los errores estándar estimados para los parámetros tenderán a ser menores que sus valores reales, conduciendo a una significación estadística erróneamente mayor. En estos casos, se utilizan mínimos cuadrados generalizados para obtener mejores estimaciones del error estándar y dar cuenta de la autocorrelación de las series residuales.
Nota: recordar que la transformación logarítmica estabiliza la varianza de la serie.
1. Modelos lineales:
Un modelo para la serie temporal {X(t), t=1,…,n} es lineal si puede expresarse como:
X(t)=alfa_0+alfa_1*U_1(t)+…+ alfa_p*U_p(t)+Z(t)
donde U_i(t) es el valor de la variable explicativa o predictor ith en el tiempo t, Z(t) es el error en el tiempo t y alfa_i son los parámetros del modelo, que pueden ser estimados por mínimos cuadrados.Ejemplo: función polinómica de orden p.
X(t)=alfa_0+alfa_1*t+alfa_2*t^2+…+alfa_p*t^p+Z(t)Nota: muchos modelos no-lineales pueden transformarse en modelos lineales. Por ejemplo, el modelo X(t)=exp(alfa_0+alfa_1*t+Z(t)) de la serie {X(t)} puede transformarse tomando logaritmos naturales para obtener un modelo lineal de la serie {Y(t)} mediante Y(t)=log(X(t)). Luego, para obtener predicciones sobre X(t), debemos utilizar la transformación inversa sobre Y(t): exp(Y(t)). (ver sección 5.10 para una discusión de este tema)
Nota: la diferenciación de la serie puede eliminar la tendencia estocástica y determinística: diff(X(t))=X(t)-X(t-1). Por ejemplo, si la tendencia subyacente es polinómica de orden m, necesitaremos una diferenciación de orden m para eliminar la tendencia.
- Mínimos Cuadrados Ordinarios (OLS)
Aplicación en R: Mínimos Cuadrados simples (OLS).
#Ajuste de un modelo lineal para datos simulados: simulamos una serie donde los errores están autocorrelacionados: X(t)=50+3*t+Z(t)
library(fma)
set.seed(1)
z <- w <- rnorm(100, sd = 20) #errores
for (t in 2:100) z[t] <- 0.8 * z[t - 1] + w[t] #errores autocorrelacionados
Time <- 1:100
x <- 50 + 3 * Time + z #serie simulada
plot(x, xlab = "time", type = "l") #gráfico de la serie
x.lm <- lm(x ~ Time) #ajuste del modelo lineal
coef(x.lm) # coeficientes estimados
sqrt(diag(vcov(x.lm))) #errores estándar (sub-estimados por la autocorrelación de los residuales)
summary(x.lm) #resumen del ajuste
acf(resid(x.lm)) #como esperábamos, la serie temporal de los residuos está autocorrelacionada, sugiriendo que los residuos siguen un proceso AR(1) (lo simulamos así).
#Ajuste de modelo lineal para la serie de temperaturas (1970-2005): ajuste de un modelo lineal y construcción de los IC al 95% para los parámetros estimados.
www <- "http://www.massey.ac.nz/~pscowper/ts/global.dat"
Global <- scan(www)
Global.ts <- ts(Global, st = c(1856, 1), end = c(2005,12), fr = 12)
temp <- window(Global.ts, start = 1970)
temp.lm <- lm(temp ~ time(temp)) #ajustamos un modelo lineal
coef(temp.lm) #parámetros estimados
confint(temp.lm) #IC al 95% para los parámetros estimados: no contienen el 0, lo que da evidencia estadística de una tendencia creciente de las temperaturas globales si la autocorrelación de los residuales es nula.
acf(resid(lm(temp ~ time(temp)))) #ACF para los residuos del modelo: la serie está positivamente autocorrelacionada en lags bajos, lo que indica una sub-estimación del error estandar e IC muy estrechos para la pendiente.
- Mínimos Cuadrados Generalizados (GLS)
Aplicación en R: Mínimos Cuadrados Generalizados (GLS).
#Ajuste de un modelo lineal por GLS para datos simulados: simulamos una serie donde los errores están autocorrelacionados: X(t)=50+3*t+Z(t)
library(nlme)
x.gls <- gls(x ~ Time, cor = corAR1(0.8)) #ajuste de modelo lineal por GLS
coef(x.gls) #parámetros estimados
sqrt(diag(vcov(x.gls))) #errores estándar para los parámetros estimados
acf(resid(x.gls))
#los errores estándar de los parámetros son mayores que los obtenidos por OLS (lm) y por tanto más adecuados, ya que tienen en cuenta la autocorrelaqción.
#Ajuste de modelo lineal por GLS para la serie de temperaturas (1970-2005): ajuste de un modelo lineal y construcción de los IC al 95% para los parámetros estimados.
#necesitamos estimar la autocorrelación en lag=1 mediante OLS para incluirlo en el GLS
acf(resid(temp.lm)) #la serie residual parece segui un proceso AR(1) con autocorrelación en lag=1 de 0.7.
temp.gls <- gls(temp ~ time(temp), cor = corAR1(0.7)) #ajuste de modelo lineal por GLS. Se utiliza la autocorrelación anterior como parámetro.
coef(temp.gls) #los parámetros estimados son diferentes que los obtenidos por OLS debido a la ponderación.
confint(temp.gls) #los IC al 95% son más anchos que los obtenidos por OLS, pero el 0 no está contenido en ellos, lo que implica que las estimaciones son estadísticamente significativas y por tanto, la tendencia es significativa.
#en conclusión, existe evidencia estadística de una tendencia en aumento de la temperatura global en el período 1970-2005, por lo que si las condiciones actuales persisten, es de esperar que las temperaturas continúen aumentando en el futuro.
2. Modelo lineal con variables estacionales
En las series temporales pueden existir efectos estacionales debido especialmente a los ciclos anuales causados directa o indirectamente por el movimiento de la Tierra alrededor del Sol.
- Variables indicadoras estacionales aditivas
Ejemplo: para una serie temporal {X(t)} mensual (t=1, corresponde al mes de Enero), podemos estimar los parámetros del siguiente modelo mediante OLS o GLS, tratando al término estacional S(t) como factor para obtener los índices estacionales.
Es decir, el modelo es X(t)=alfa_1+S(t)+Z(t), donde S(t) es un factor con parámetros beta_1 (para t=1,13,…),…,beta_12 (para t=12,24,…).
Aplicación en R: Serie de Temperaturas (1970-2005)
### Regresión lineal con variables indicadoras estacionales ###
#estimamos los parámetros de la tendencia lineal con índices estacionales aditivos para la serie temporal de 1970-2005.
Seas <- cycle(temp) #extrae los índices estacionales
Time <- time(temp)
#estimación con OLS
temp.lm <- lm(temp ~ 0 + Time + factor(Seas)) #estimamos los parámetros para el modelo con indicadores estacionales mediante OLS tratando al término estacional S(t) como factor y poniendo "0" dentro de la fórmula para asegurar que el modelo no tiene intecepto.
#si incluyéramos el intercepto en el modelo, uno de los términos estacional sería eliminado y en la salida veeríamos una estimación del intercepto. Los resultados del ajuste serían equivalentes.
coef(temp.lm) #parámetros estimados
#predicción para los 2 años siguientes
new.t <- seq(2006, len = 2 * 12, by = 1/12)
alpha <- coef(temp.lm)[1]
beta <- rep(coef(temp.lm)[2:13], 2)
(alpha * new.t + beta)[1:4]
#alternativamente podemos realizar la predicción con las etiquetas estacionales correctas
new.dat <- data.frame(Time = new.t, Seas = rep(1:12, 2))
predict(temp.lm, new.dat)[1:24]
- Variables estacionales armónicas
Generalmente, se utilizan las funciones seno y coseno para construir la variación estacional del modelo estacional.
Ejemplo: Una ola seno con frecuencia f (ciclos por intervalos muestreado), de amplitud A y cambio de fase phi, se expresa como:
A*sen(2*phi*f*t+phi)=alfa_s*sen(2*phi*f*t)+alfa_c*cos(2*phi*f*t)
Por tanto, para una serie temporal {X(t)} con s períodos estacionales existen [s/2] ciclos posibles, el modelo de variables armónicas es:X(t)=m(t)+sum(s_i*sin(2*^phi*i*t/s)+c_i*cos(2*phi*i*t/s))+Z(t)
donde la suma va desde i=1 a i=[s/2], m(t) es la tendencia que incluye un parámetro para el término constante y s_i y c_i son parámetros desconocidos. Cuando s es un número par, el seno se vuelve nulo. Entonces, este modelo armónico con el término constante tendrá el mismo número de parámetros máximo que el modelo con variables indicadoras estacionales.Aplicación en R: simularemos el modelo anterior de variables estacionales armónicas. Por ejemplo, para: X(t)= 0.1 + 0.005*t + 0.001*t^2 + sin(2*phi*t/12)+0.2*sin(4*phi*t/12) + 0.1*sin(8*phi*t/12) + 0.1*cos(8*phi*t/12) + W(t). con W(t) un ruido blanco gaussiano con desviación estándar .5.
### Regresión lineal con variables estacionales armónicas ###
#Ajustamos un modelo simulado: X(t)= 0.1 + 0.005*t + 0.001*t^2 + sin(2*phi*t/12)+0.2*sin(4*phi*t/12) + 0.1*sin(8*phi*t/12) + 0.1*cos(8*phi*t/12) + W(t). con W(t) un ruido blanco gaussiano con desviación estándar .5.
set.seed(1)
TIME <- 1:(10 * 12)
w <- rnorm(10 * 12, sd = 0.5)
Trend <- 0.1 + 0.005 * TIME + 0.001 * TIME^2
Seasonal <- sin(2*pi*TIME/12) + 0.2*sin(2*pi*2*TIME/12) +0.1*sin(2*pi*4*TIME/12) + 0.1*cos(2*pi*4*TIME/12)
x <- Trend + Seasonal + w
plot(x, type = "l")
#ponemos en matrices las variabels armónicas
SIN <- COS <- matrix(nr = length(TIME), nc = 6)
for (i in 1:6)
{
COS[, i] <- cos(2 * pi * i * TIME/12)
SIN[, i] <- sin(2 * pi * i * TIME/12)
}
#en general, no conocemos el orden de la tendencia armónica y polinómica, y se vuelve subjetivo decidir qué variables son significativas. El t-ratio de magnitud 2 es una forma común de elección y corresponde a un nivel de significación aprox del 5%.
x.lm1 <- lm(x ~ TIME + I(TIME^2) + COS[, 1] + SIN[, 1] +COS[, 2] + SIN[, 2] + COS[, 3] + SIN[, 3] + COS[, 4] +SIN[, 4] + COS[, 5] + SIN[, 5] + COS[, 6] + SIN[, 6])
coef(x.lm1)/sqrt(diag(vcov(x.lm1))) #obtenemos el t-ratio dividiendo los coeficientes por los errores estándares de las estimaciones.
#se observan 3 coeficientes significativos, que serán los utilizados en el siguiente modelo
x.lm2 <- lm(x ~ I(TIME^2) + SIN[, 1] + SIN[, 2])
coef(x.lm2)/sqrt(diag(vcov(x.lm2))) #obtenemos el t-ratio dividiendo los coeficientes por los errores estándares de las estimaciones.
#se observa que todos los coeficientes son significativos.
coef(x.lm2) #coeficientes estimados para el modelo mejor ajustado
#los coeficientes dan el siguiente modelo para las predicciones en el tiempo t: xt = 0.280 + 0.00104t2 + 0.900 sin(2 t/12) + 0.199 sin(4 t/12)
#comparamos los 2 modelos ajustados mediante el criterio de AIC y el modelo generador de las simulaciones
AIC(lm(x ~ TIME +I(TIME^2) +SIN[,1] +SIN[,2] +SIN[,4] +COS[,4])) #AIC del modelo subyacente
AIC(x.lm1)
AIC(x.lm2)
#observamos que el último modelo tiene el menor AIC y por tanto da el mejor ajuste a los datos.
#debido a la variación de muestreo el modelo mejor ajustado no es idéntico al modelo utilizado para simular los datos.
#podemos automatizar la selección del mejor modelo por AIC utilizando la función "step"
step(x.lm1)
#el mejor ajuste puede elegirse según el modelo que conduzca a la menor desviación estandar estimada para los errores, según los grados de libertad considerados.
#Ajustamos un modelo armónico para las series de temperatura, con una tendencia cuadrática.
#como la variable está en anos, no necesitamos el divisor "12" para crear las variables armónicas.
#estandarizamos la variable TIME luego de calcular las predictoras COS y SIn para reducir el cálculo del OLS.
SIN <- COS <- matrix(nr = length(temp), nc = 6)
for (i in 1:6)
{
COS[, i] <- cos(2 * pi * i * time(temp))
SIN[, i] <- sin(2 * pi * i * time(temp))
}
TIME <- (time(temp) - mean(time(temp)))/sd(time(temp))
mean(time(temp))
sd(time(temp))
temp.lm1 <- lm(temp ~ TIME + I(TIME^2)
+COS[,1] + SIN[,1] + COS[,2] + SIN[,2]
+COS[,3] + SIN[,3] + COS[,4] + SIN[,4]
+COS[,5] + SIN[,5] + COS[,6] + SIN[,6])
coef(temp.lm1)/sqrt(diag(vcov(temp.lm1)))
temp.lm2 <- lm(temp ~ TIME + SIN[, 1] + SIN[, 2])
coef(temp.lm2)
AIC(temp.lm)
AIC(temp.lm1)
AIC(temp.lm2)
#para chequear el modelo ajustado creamos el gráfico del tiempo (para detectar patrones en la serie) y el correlograma (para detectar autocorrelación) de los residuales.
plot(time(temp), resid(temp.lm2), type = "l")
abline(0, 0, col = "red")
acf(resid(temp.lm2))
pacf(resid(temp.lm2))
#no observamos curvas en la serie, lo que implica que la línea recta es una descripción adecuada para la tendencia. La persistencia de la tendencia en estar encima y debajo del eje-x indica que la serie está positivamente autocorrelacionada.
#verificamos la afirmación anterior con el correlograma de los residuales, que indica una autocorrelación positiva en lag=1-10. Además, es similar al esperado para un proceso AR(p).
#verificamos el modelo propuesto con el PACF, donde solo los lags=1 y 2 son estadísticamente significativos (indica AR(2)).
#ajustamos un modelo AR(2)
res.ar <- ar(resid(temp.lm2), method = "mle")
res.ar$ar
sd(res.ar$res[-(1:2)])
acf(res.ar$res[-(1:2)]) #los residuales son aproximadamente ruido blanco
#la forma final del modelo ajustado para la ST de temperatura mensual es: x(t) = 0.175 +0.184*(t - 1988)/10.4+ 0.0204*sin(2*phi*t) + 0.0162*sin(4*phi*t) + z(t), y {Z(t)} sigue un proceso AR(2) dado por z(t) = 0.494*z(t-1) + 0.307*z(t-2) + w(t) donde {W(t)} es un ruido blanco con media cero y desviación estándar .0837.
- Transformación logarítmica
Para series con valores negativos, podemos sumarle una constante a cada término de la serie, de tal manera que si {X(t)} es la serie que contiene algunos valores negativos, si le agregamos una constante c_0 tal que c_0>max{-X(t)} y tomamos logaritmos, obtenemos una serie transformada {log(c_0+X(t))} que está definida para todo t. Entonces, podemos ajustar un modelo lineal (e.g. tendencia lineal): x(t)=-c_0+exp(alfa_0+alfa_1*t+z(t)), donde alfa_0 y alfa-1 son los parámetros del modelo y {z(t)} es la serie residual (que puede estar autocorrelacionada). El problema es que el parámetro c_0 también debe ser estimado por el modelo, cuando la selección de c_0 suele hacerse tal que satisfaga la condición c_0>max(-X(t)). Si no existe evidencia para términos residuales multiplicativos, la constante c_0 será estimada con los demás parámetros del modelo, utilizando mínimos cuadrados no-lineales, ajustando el siguiente modelo: x(t)=-c_0+exp(alfa_0+alfa_1*t+z(t)).
3. Modelos no-lineales (e.g. modelo de Bass):
Para algunas series (e.g. cuando conocemos que el proceso subyacente es no-lineal) es más apropiado ajustar un modelo no-lineal directamente más que tomar logaritmos o utilizar aproximaciones polinómicas lineales.
Aplicación en R: simularemos una serie de datos no-lineal con c_0=0 y comparamos los parámetros estimados utilizando la función “nls” con aquellos parámetros conocidos de la función subyacente.
#### Regresión con modelos no-lineales ####
#simulamos una serie no lineal con residuales AR(1)
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) #la serie tiene una tendencia exponencial creciente pero además contiene valores egativos, por loque no podemos utilizar directamente la transformación logarítmica, sino que deberemos utilizar el modelo no-lineal.
x.nls <- nls(x ~ exp(alp0 + alp1 * Time), start = list(alp0 = 0.1,alp1 = 0.5)) #ajustamos un modelo no-lineal especificando la fórmula con los parámetros y sus valores de comienzo.
summary(x.nls)$parameters #los valores estimados para alfa_0 y alfa_1 son cercanos a los valores subyacentes que se utilizaron para simular los datos, pero los errores estandar para las estimaciones aprece estar sub-estimados porque existe autocorrelación de los residuales.
acf(resid(x.nls))
Comentarios
Publicar un comentario