Tema 7: Modelos de regresión dinámica

shuwei.chou@ucr.ac.cr

Escuela de Estadística, Universidad de Costa Rica

Paquetes de R

Para este tema, se necesita cargar estos paquetes:

library(ggplot2)
library(forecast)
library(fpp2)
library(astsa)
library(sarima)

Introducción

Contenido

Introducción
Regresión lineal con error de tipo ARIMA
Serie estacionaria por tendencia y por diferencia.
Análisis de intervención

Introducción

En el Tema 4, hemos visto regresión con series de tiempo:

\[Y_t=\beta_0+\beta_1 X_{t,1}+\beta_2 X_{t,2}+...+\beta_p X_{t,p}+\epsilon_t, ~~t=1,...,T,\] en donde \(\left\lbrace \epsilon_t \right\rbrace\) es ruido blanco Gaussiano.

Los modelos de regresión dinámicos son modelos de regresión en los que el término de error sigue un modelo ARIMA, y la variable independiente \(X_t\) influye a \(Y_t\) no solo instantáneamente en el tiempo \(t\) sino a través de varios periodos.
El caso univariado, \[Y_t=\beta_0+\beta_1 X_{t}+\beta_2 X_{t-1}+...+\beta_p X_{t-p}+\eta_t\] \[\left\lbrace \eta_t \right\rbrace \sim ARIMA(p,d,q)\]
Veremos algunas extensiones del modelo de regresión con error tipo ARIMA antes de considerar la variable independiente rezagada como predictora.

Regresión lineal con error de tipo ARIMA

Contenido

Introducción
Regresión lineal con error de tipo ARIMA
Serie estacionaria por tendencia y por diferencia.
Análisis de intervención

Regresión lineal con error de tipo ARIMA

Primeramente, vamos a considerar el siguiente modelo: \[Y_t=\beta_0+\beta_1 X_{t,1}+\beta_2 X_{t,2}+...+\beta_p X_{t,p}+\eta_t, ~~t=1,...,T,\] en donde \(\left\lbrace \eta_t \right\rbrace \sim ARIMA(p,d,q)\), i.e.

\[\phi(B) (1-B)^d \eta_t=\theta(B) \epsilon_t,\] con ruido blanco \(\left\lbrace \epsilon_t \right\rbrace\).

La estimación de los parámetros del modelo minimizando la suma de cuadrados \(\eta_t\) en vez de \(\eta_t\) ignorando su estructura dependiente del tiempo es incorrecta.
La estimación de los coeficientes \(\beta_i, i=1,...,p\) no satisfacen las propiedades óptimas (insesgamiento, variancia mínima, etc.). Por lo tanto, toda la teoría de pruebas de hipótesis de regresión no funciona.

Notas:

Verifique que la variable \(Y_t\) y las \(X_{t1},X_{t2},...,X_{tp}\) sean estacionarias, porque si se estiman los coeficientes con variables no estacionarias, los estimadores no son consistentes.
Además, podría presentar problema de la regresión espuria.
Realice diferencias a las variables no estacionarias.
El modelo después de aplicar la diferencia es llamado modelo de regresión en diferencias.
Se puede demostrar que el modelo de regresión con error de tipo ARIMA es equivalente a un modelo de regresión en diferencias con error tipo ARMA.

Por ejemplo, considere \[Y_t=\beta_0+\beta_1 X_{t,1}+\beta_2 X_{t,2}+...+\beta_p X_{t,p}+\eta_t, ~~t=1,...,T,\] \[(1-\phi_1B) (1-B) \eta_t=(1-\theta_1 B) \epsilon_t,\]
Después de aplicar las diferencias \[Y'_t=Y_{t}-Y_{t-1}, ~~~~ X'_{t,k}=X_{t,k}-X_{t-1,k},~k=1,...,p,~~\text{y}~~ \eta'_t=\eta_t-\eta_{t-1}\]

obtenemos

\[Y'_t=\beta_1 X'_{t,1}+\beta_2 X'_{t,2}+...+\beta_p X'_{t,p}+\eta'_t, ~~t=1,...,T,\] \[(1-\phi_1B) \eta'_t=(1-\theta_1 B) \epsilon_t,\]

En R, si queremos ajustar un modelo de regresión con error tipo ARIMA(1,1,0), i.e. \[Y_t=\beta_0+\beta_1 X_{t,1}+\eta_t\] \[(1-\phi_1B) (1-B) \eta_t=\epsilon_t,\] con el siguiente comando:

mod <- Arima(y, xreg=x, order=c(1,1,0))

El programa considera el modelo en diferencias con error tipo AR(1):

\[Y'_t=\beta_1 X'_{t1}+\eta_t\] \[(1-\phi_1B) \eta_t=\epsilon_t,\]

Ejemplo tomado de Hyndman (2018): pronóstico del cambio de gasto basado en el ingreso personal (serie trimestral) de 01-1970 a 03-2016.

autoplot(uschange[,1:2], facets=TRUE) +
  xlab("Year") + ylab("") +
  ggtitle("Quarterly changes in US consumption
    and personal income")

(mod <- auto.arima(uschange[,"Consumption"],
                   xreg=uschange[,"Income"]))

Series: uschange[, "Consumption"] 
Regression with ARIMA(1,0,2) errors 

Coefficients:
         ar1      ma1     ma2  intercept    xreg
      0.6922  -0.5758  0.1984     0.5990  0.2028
s.e.  0.1159   0.1301  0.0756     0.0884  0.0461

sigma^2 = 0.3219:  log likelihood = -156.95
AIC=325.91   AICc=326.37   BIC=345.29

El modelo final estimado es:

\[Y_t=0.599+ 0.203 X_{t1}+\eta_t\] \[\eta_t=0.692 \eta_{t-1}+\epsilon_t-0.576 \epsilon_{t-1}+ 0.198\epsilon_{t-2},\] \[\epsilon_{t} \sim N(0,0.322)\]

errores <- cbind("Regression Errors" = 
                   residuals(mod, type="regression"),
                 "ARIMA errors" = 
                   residuals(mod, type="innovation")) 

head(errores)

        Regression Errors ARIMA errors
1970 Q1       -0.18024704  -0.16714211
1970 Q2       -0.37577719  -0.31981056
1970 Q3       -0.03728173   0.07199692
1970 Q4       -0.82151101  -0.69355381
1971 Q1        0.89529777   1.05009349
1971 Q2        0.01940569   0.14169806

errores%>%
  autoplot(facets=TRUE)

Residuales del modelo de regresión con errores independientes.

acf2(residuals(mod, type="regression"))

Residuales del modelo dinámico con error de estructura ARMA.

acf2(residuals(mod, type="innovation"))

checkresiduals(mod)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(1,0,2) errors
Q* = 5.8916, df = 5, p-value = 0.3169

Model df: 3.   Total lags used: 8

checkresiduals(mod, lag = 16)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(1,0,2) errors
Q* = 14.068, df = 13, p-value = 0.3691

Model df: 3.   Total lags used: 16

Pronóstico

Se necesitan valores futuros de la covariable.
Una posibilidad es suponer que toma el promedio histórico.

fcast <- forecast(mod, xreg=rep(mean(uschange[,2]),8))
autoplot(fcast) + xlab("Year") +
  ylab("Percentage change")

Serie estacionaria por tendencia y por diferencia.

Contenido

Introducción
Regresión lineal con error de tipo ARIMA
Serie estacionaria por tendencia y por diferencia.
Análisis de intervención

Serie estacionaria por tendencia y por diferencia.

Suponga que una serie temporal \(\left\lbrace Y_t \right\rbrace\) es una realización de una tendencia determinística y un componente estocástico:

\[ Y_t=TD_t+ \eta_t, \] donde \(TD_t=\beta_0+\beta_1 t\) y \(\eta_t \sim ARIMA(p,d,q)\).

Caso 1: si \(d=0\), \(\left\lbrace Y_t \right\rbrace\) es estacionaria alrededor de una tendencia determinística. Por lo tanto, se puede eliminar la tendencia de la serie original y ajustar un modelo ARMA a los residuales.

Caso 2: si \(d>0\), \(\left\lbrace Y_t \right\rbrace\) es estacionaria por diferencia. Por lo tanto, se puede realizar una diferencia para obtener una serie estacionaria. Caso más común es cuando \(d=1\).

Ejemplo de estos dos tipos de estacionariedad:

Tendencia determinística: \[Y_t=Y_{t-1}+\mu=Y_0+\mu t\] Tendencia estocástica (acumulación de choques aleatorias): \[Y_t=Y_{t-1}+\epsilon_t=Y_0+\sum_{s=1}^t \epsilon_s\] donde \(\mu\) es una constante y \(\epsilon_t\) es ruido blanco.

En síntesis, una serie temporal \(\left\lbrace Y_t \right\rbrace\) está compuesto por una tendencia determinística y un componente estocástico que es modelado por \(ARIMA(p,d,q)\).
Se puede descomponer \(\eta_t\) en dos componentes: tendencia estocástica (choques aleatorios) y el componente aleatorio “estacionario”.
Entonces, \(\left\lbrace Y_t \right\rbrace\) se puede descomponer en tres componentes:

tendencia determinística,
tendencia estocástica, y
el componente “aleatorio”.

Un modelo estacionario por tendencia, no tiene la tendencia estocástica, y el componente aleatorio es \(ARMA(p,q)\).
En el caso de un modelo estacionario por diferencia, el polinomio autoregresivo del componente \(\eta_t\) tiene al menos una raíz unitaria.

Tendencia determinística: \(Y_t= 0 + 0.5 t + \epsilon_t\)

Tendencia determinística con tendencia estocástica: \(Y_t= 0 + 0.5 t + \sum_{s=1}^t \epsilon_s\)

Tendencia estocástica: \(Y_t = \sum_{s=1}^t \epsilon_s\)

Tendencia determinística y estocástica

Devolviendo al modelo regresión, en la práctica se puede modelar una tendencia lineal usando: \[Y_t=\beta_0+\beta_1 t + \eta_t\]
1. \(\eta_t \sim ARMA(p,q)\), o
2. \(\eta_t \sim ARIMA(p,1,q)\).
En el caso 2, se puede simplificar el modelo en: \[Y_t=Y_{t-1}+\beta_1+ \eta'_t.\] Este modelo es similar a un modelo de caminata aleatoria pero con un desvío \(\beta_1\) y el error es ARMA.

autoplot(austa) + xlab("Year") +
  ylab("millions of people") +
  ggtitle("Total annual international visitors to Australia")

Ajuste con tendencia determinística

trend <- seq_along(austa)
(fit1 <- auto.arima(austa, d=0, xreg=trend))

Series: austa 
Regression with ARIMA(2,0,0) errors 

Coefficients:
         ar1      ar2  intercept    xreg
      1.1127  -0.3805     0.4156  0.1710
s.e.  0.1600   0.1585     0.1897  0.0088

sigma^2 = 0.02979:  log likelihood = 13.6
AIC=-17.2   AICc=-15.2   BIC=-9.28

\[Y_t=0.416+0.171t+\eta_t\] \[\eta_t=1.113\eta_{t-1}-0.380 \eta_{t-2}+\epsilon_t\] \[\epsilon_t \overset{iid}{\sim} N(0,0.03)\]

Ajuste con tendencia estocástica

(fit2 <- auto.arima(austa, d=1))

Series: austa 
ARIMA(0,1,1) with drift 

Coefficients:
         ma1   drift
      0.3006  0.1735
s.e.  0.1647  0.0390

sigma^2 = 0.03376:  log likelihood = 10.62
AIC=-15.24   AICc=-14.46   BIC=-10.57

\[Y_t-Y_{t-1}=0.173+\eta'_t,\] o de otra forma,

\[Y_t=Y_0+0.173t+\eta_t\] \[\eta_t=\eta_{t-1}+0.301\epsilon_{t-1}+\epsilon_t\] \[\epsilon_t \overset{iid}{\sim} N(0,0.034)\]

fc1 <- forecast(fit1,
  xreg = length(austa) + 1:10)
fc2 <- forecast(fit2, h=10)
autoplot(austa) +
  autolayer(fc2, series="Stochastic trend") +
  autolayer(fc1, series="Deterministic trend") +
  ggtitle("Forecasts from trend models") +
  xlab("Year") + ylab("Visitors to Australia (millions)") +
  guides(colour=guide_legend(title="Forecast"))

Análisis de intervención

Contenido

Introducción
Regresión lineal con error de tipo ARIMA
Serie estacionaria por tendencia y por diferencia.
Análisis de intervención

Análisis de intervención

Una serie puede experimentar cambios de comportamiento en el tiempo cuando ocurren fenómenos como cambios en las políticas públicas, desastres naturales, crisis económicas serias, aumentos sustanciales en gastos de publicidad, etc.
La suposición del mismo comportamiento a lo largo de tiempo puede ser no realista.
El análisis de intervención introducida por Box & Tiao en 1975 toma en cuenta este tipo de eventos.
Ilustramos el análisis con el siguiente ejemplo (Wichern & Jones, 1977).

Series semanales de cuotas de mercado dentífrico de las marcas Colgate y Crest en los Estados Unidos de 1 de enero de 1958 al abril de 1963.
Al inicio, Colgate aventajaba a la marca Crest en el mercado.
El 1 de agosto de 1960 ocurrió el cambio de comportamiento: la Asociación Dental Americana dio un respaldo enorme a la marca Crest al hacer público que era una pasta dental eficaz para prevenir las caries dentales.
Procter y Gamble, la compañía que producía la marca Crest, aprovechó y divulgó intensamente durante dos semanas el enuncio.

table<-read.table("crestcolgate.dat")
colnames(table)<-c( " CRESTMS", "COLGTEMS" ,"CRESTPR", "COLGTEPR")

crest.colgate<-ts(table[,1:2])
autoplot(crest.colgate)

Pasos a seguir:

Elaborar un modelo ARIMA antes de la intervención, i.e. durante las primeras 134 semanas.
Ampliar el modelo agregando variables indicadoras para registrar la intervención.
Re-estimar el modelo con la serie completa con las variables indicadoras.

crest<-ts(table[,1])
crest1<-crest.colgate[1:134]
crest2<-crest.colgate[135:276]

ts.plot(crest)
points(135:276,crest2,col=2,type="l")
legend("topleft",c("antes","después"),lwd=2,col=1:2)

Crest1: antes de la intervención

acf2(crest1,main="antes")

Serie diferenciada (antes de la intervención)

dif.crest1<-diff(crest1)
ts.plot(dif.crest1)

acf2(dif.crest1)

Vamos a ajustar un ARIMA(0,1,1): \[(1-B)Z_t=(1-\theta_1B)a_t\]

moda <- Arima(crest1, order=c(0,1,1))
summary(moda)

Series: crest1 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.6918
s.e.   0.0644

sigma^2 = 0.001248:  log likelihood = 256.11
AIC=-508.22   AICc=-508.12   BIC=-502.44

Training set error measures:
                       ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.0006404213 0.03505771 0.02747198 -5.725771 22.02557 0.8434379
                      ACF1
Training set -0.0009184472

El modelo se puede despejar \(Z_t\): \[Z_t=\frac{(1-\theta_1B)}{(1-B)}a_t\]
Creación de variables indicadoras para indicar la intervención: \[I_{1t}=\left\lbrace \begin{aligned} 0, & & \text{si}~~ t<135 \\ 1, & & \text{si}~~ t \geq135, \end{aligned} \right., ~~ \text{y}~~~I_{2t}=\left\lbrace \begin{aligned} 0, & & \text{si}~~ t<136 \\ 1, & & \text{si}~~ t \geq136, \end{aligned} \right.\]

I1 <- c(rep(0,134),rep(1,142)) #intervención durante la semana 135
I2 <- c(rep(0,135),rep(1,141)) #intervención durante la semana 136
X<-cbind(I1,I2)

X[130:140,]

      I1 I2
 [1,]  0  0
 [2,]  0  0
 [3,]  0  0
 [4,]  0  0
 [5,]  0  0
 [6,]  1  0
 [7,]  1  1
 [8,]  1  1
 [9,]  1  1
[10,]  1  1
[11,]  1  1

El modelo ARIMA con intervención: \[Z_t=c_1 I_{1t} +c_2 I_{2t}+ \frac{(1-\theta_1B)}{(1-B)}a_t\]
Multiplicando por ambos lados \((1-B)\): \[(1-B) Z_t=c_1 (1-B) I_{1t} +c_2 (1-B) I_{2t}+ (1-\theta_1B)a_t\] Note que en realidad…

X_diff <- cbind(diff(I1),diff(I2))
X_diff[130:140,]

      [,1] [,2]
 [1,]    0    0
 [2,]    0    0
 [3,]    0    0
 [4,]    0    0
 [5,]    1    0
 [6,]    0    1
 [7,]    0    0
 [8,]    0    0
 [9,]    0    0
[10,]    0    0
[11,]    0    0

modb <- Arima(crest, xreg=X, order=c(0,1,1))
summary(modb)

Series: crest 
Regression with ARIMA(0,1,1) errors 

Coefficients:
          ma1      I1      I2
      -0.7782  0.0654  0.1119
s.e.   0.0437  0.0434  0.0434

sigma^2 = 0.001902:  log likelihood = 472.22
AIC=-936.44   AICc=-936.29   BIC=-921.98

Training set error measures:
                     ME       RMSE        MAE      MPE     MAPE      MASE
Training set 0.00186971 0.04329952 0.03363337 -3.49533 16.66087 0.7903928
                    ACF1
Training set -0.01155037

checkresiduals(modb,lag=20)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(0,1,1) errors
Q* = 15.279, df = 19, p-value = 0.7047

Model df: 1.   Total lags used: 20

Otros tipos de intervención

Capítulo 11 del Cryer & Chan (2008) Time Series Analysis with Applications in R.

Tema 7: Modelos de regresión dinámica

Paquetes de R

Introducción

Introducción

Regresión lineal con error de tipo ARIMA

Regresión lineal con error de tipo ARIMA

Pronóstico

Serie estacionaria por tendencia y por diferencia.

Serie estacionaria por tendencia y por diferencia.

Tendencia determinística y estocástica

Ajuste con tendencia determinística

Ajuste con tendencia estocástica

Análisis de intervención

Análisis de intervención

Otros tipos de intervención

Próximo tema

Tema 8: ARCH/GARCH