5.3 Modelo lineal con series temporales
En las secciones anteriores se ha visto como cargar y crear un objeto de clase zooreg para tratar datos de series temporales. A continuación, se presentan algunas funciones específicas para la estimación del modelo lineal general en incluyendo elementos habituales para esta clase de datos como retardos y diferencias. Tomando el conjunto de datos sobre consumo de tabaco y publicidad, cuya descripción se puede encontrar al final, se cargan los datos y se crea un objeto zooreg excluyendo la primera columna cargada (que contiene el año):
> library(zoo)
> tabd <- read.csv("http://grserrano.es/datos/datos_tabaco.csv")
> tabdz <- zooreg(tabd[, -1], start=1930, frequency=1)
La falta de estacionariedad (media no constante) de las variables del ejemplo es evidente a la vista del gráfico en logaritmos que se presentó en la entrada sobre gráficos de series temporales, por lo que especificaremos los modelos en tasas de variación logarítmica. Para justificar formalmente tomar diferencias, se pueden utilizar los contrastes de Dickey-Fuller (adf.test(), que se ilustra un poco más adelante) o el de Phillips-Perron (pp.test()) aplicados a las series en logaritmos.
Para la estimación, el paquete dynlm estima modelos lineales por mínimos cuadrados, pero facilita la especificación de regresiones con retardos y diferencias en las variables. No insistiré aquí en la importancia de no estimar regresiones con variables que comparten tendencias, puesto que casi siempre serán regresiones espurias.
En el siguiente fragmento se estima un modelo en tasas de variación (diferencias de logartimos, usando la función d()) explicando el consumo de tabaco en función del gasto en publicidad utilizando datos desde 1940 (argumento data con la función window()) que se muestran en el Gráfico 1.
Gráico 1: Evolución del consumo de tabaco y el gasto en publicidad en tasas logarítmicas de variación interanual.
> library(dynlm)
> tabdz1 <- diff(log(window(tabdz[, c("consumo", "publicidad")], start=1940)))
> tmp <- data.frame(date=as.Date(paste(1941:1978, "-6-30", sep="")),
+ consumo=tabdz1[, "consumo"],
+ publicidad=tabdz1[, "publicidad"])
> tmpl <- melt(tmp, id.vars="date")
> qplot(date, value*100, data=tmpl, geom=c("point", "line"),
+ ylab="T.V. Interanual (%)", xlab="Año", colour=variable)
> tabm1 <- dynlm(d(log(consumo))~d(log(publicidad)),
+ data=window(tabdz, start=1940))
> summary(tabm1)
Time series regression with "zooreg" data:
Start = 1941, End = 1978
Call:
dynlm(formula = d(log(consumo)) ~ d(log(publicidad)), data = window(tabdz,
start = 1940))
Residuals:
Min 1Q Median 3Q Max
-0.077950 -0.032970 -0.004458 0.012886 0.098907
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.026705 0.007478 3.571 0.00103 **
d(log(publicidad)) -0.262784 0.087212 -3.013 0.00471 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04285 on 36 degrees of freedom
Multiple R-squared: 0.2014, Adjusted R-squared: 0.1792
F-statistic: 9.079 on 1 and 36 DF, p-value: 0.004712
Como se aprecia en la estimación, el consumo está inversamente relacionado con el gasto publicitario, resultado anómalo sólo explicable por la progresiva concienciación de los perjuicios del tabaco para la
salud.
Para introducir un retardo de la tasa de variación del consumo de tabaco como variable explicativa se usa la función L() (dentro de la fórmula de dynlm()):
> tabm2 <- dynlm(d(log(consumo))~L(d(log(consumo))) + d(log(publicidad)),
+ data=window(tabdz, start=1940))
> summary(tabm2)
Time series regression with "zooreg" data:
Start = 1942, End = 1978
Call:
dynlm(formula = d(log(consumo)) ~ L(d(log(consumo))) + d(log(publicidad)),
data = window(tabdz, start = 1940))
Residuals:
Min 1Q Median 3Q Max
-0.07860 -0.02460 -0.00015 0.02212 0.08261
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.011739 0.007496 1.566 0.12659
L(d(log(consumo))) 0.430304 0.138768 3.101 0.00386 **
d(log(publicidad)) -0.136513 0.081429 -1.676 0.10282
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03638 on 34 degrees of freedom
Multiple R-squared: 0.3658, Adjusted R-squared: 0.3285
F-statistic: 9.806 on 2 and 34 DF, p-value: 0.0004341
Ahora es la tasa de variación del consumo en el año anterior lo que mejor explica la variación en el presente y el gasto publicitario ha perdido significación. No obstante, estos resultados deben tomarse con precaución hasta realizar los diagnósticos del modelo, ya que es frecuente la presencia de autocorrelación, lo que puede hacer que la estimación de las varianzas de los coeficientes resulte inconsistente (Sección 5.4).
El gráfico de residuos (Gráfico 2(a)), como siempre, es una primera herramienta para el diagnóstico de un modelo estimado. Primero hay que preparar un data.frame con las fechas y los residuos tipificados para usar ggplot2:
> tmp <- data.frame(date=as.Date(paste(1942:1978, "-6-30", sep="")),
+ resid=scale(residuals(tabm2)),
+ precio=subset(tabd, subset=(year>=1942), select="precio"))
> g6 <- qplot(date, resid, data=tmp, color=precio, geom=c("point", "line"),
+ ylab="Residuos tipificados", xlab="Año")
A la vista de las funciones de autocorrelación simple y parcial del Gráfico 2(b) (calculadas con acf() y pacf()) no se sospecha la presencia de autocorrelación en las perturbaciones del modelo, aunque sí se observa cierta variabilidad en la dispersión en el gráfico de residuos. Veremos contrastes formales en la Sección 5.4. Sobre el código del ejemplo, en la práctica no hay que programar el gráfico de las funciones de autocorrelación cada vez, sino que habría que empaquetarlo en una función que estaría siempre disponible (véase el Apéndice C).
Gráfico de residuos de tabm2
Funciones de autocorrelación simple y parcial
Para especificar más retardos de una variable se puede usar el segundo argumento de la función de retardo L() dentro de dynlm(). A continuación se estima un modelo autorregresivo de orden 2 para el consumo de tabaco:
> tabm2.ar <- dynlm(d(log(consumo))~L(d(log(consumo)), 1:2),
+ data=window(tabdz, start=1940))
> summary(tabm2.ar)
Time series regression with "zooreg" data:
Start = 1943, End = 1978
Call:
dynlm(formula = d(log(consumo)) ~ L(d(log(consumo)), 1:2), data = window(tabdz,
start = 1940))
Residuals:
Min 1Q Median 3Q Max
-0.067064 -0.018912 0.002138 0.019484 0.089672
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.003018 0.006310 0.478 0.6356
L(d(log(consumo)), 1:2)1 0.275740 0.156556 1.761 0.0875 .
L(d(log(consumo)), 1:2)2 0.213038 0.145727 1.462 0.1532
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03445 on 33 degrees of freedom
Multiple R-squared: 0.2615, Adjusted R-squared: 0.2167
F-statistic: 5.843 on 2 and 33 DF, p-value: 0.006726
Datos de consumo de tabaco y publicidad
El archivo datos_tabaco.csv contiene datos anuales desde 1930 a 1978 para EE.UU. de las siguientes variables:
-
consumo de cigarrillos per cápita,
-
renta per cápita,
-
precio medio de la cajetilla de tabaco en términos reales (es decir, deflactado por un índice de precios) en centavos,
- un indicador del gasto en
publicidad del sector tabaquero (una variable stock que se construye a partir del gasto corriente en publicidad de años anteriores, aplicando una cierta tasa de depreciación).
Además, se conoce la siguiente información:
- En 1953, la American Cancer Society y el British Medical Research Council emitieron sendos informes que mostraban cómo las tasas de mortalidad eran más elevadas entre fumadores.
- En 1964, el US General Surgeon Report concluía que fumar cigarrillos estaba ligado al cáncer de laringe.
- En 1970 se aprobó la legislación que obligaba a las compañías a mostrar en las cajetillas mensajes de advertencia sobre los peligros del tabaco.
, se trata de un caso particular. En general, la varianza de los estimadores por MCO es:
. Sin embargo, existen situaciones en las que el supuesto ideal de perturbaciones esféricas no se cumple.
por mínimos cuadrados produce un estimador inconsistente, que invalida los contrastes realizados sobre los coeficientes estimados.
son los valores ajustados del modelo en y el orden p es decisión del analista, y suele ser un valor pequeño como 3. La hipótesis nula del contraste RESET es H0:δ2=… = δp=0 y para contrastarla se puede utilizar un estadístico F y puede que sea necesario utilizar la estimación de la matriz de covarianzas consistente a la heterocedasticidad. Existe también la versión de regresores del contraste RESET, donde en la regresión auxiliar se incluyen potencias y productos cruzados de las variables explicativas en lugar de los valores ajustados del modelo de partida.
