class: center, middle, inverse, title-slide .title[ # Tema 7: Modelos de regresión dinámica ] .subtitle[ ## Curso: Análisis de series temporales ] .author[ ### Prof. Shu Wei Chou Chen ] .institute[ ### Escuela de Estadística, UCR. ] --- # Contenido 1. Introducción 2. Regresión lineal con error de tipo ARIMA 2. 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 - 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\)`. --- # Regresión lineal con error de tipo ARIMA - 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. --- # Regresión lineal con error de tipo ARIMA - 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,$$` --- # Regresión lineal con error de tipo ARIMA - 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: ```r 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,$$` --- # Regresión lineal con error de tipo ARIMA - Ejemplo tomado de Hyndman (2018): pronóstico del cambio de gasto basado en el ingreso personal (serie trimestral) de 01-1970 a 03-2016. <img src="presentacion_files/figure-html/unnamed-chunk-3-1.png" width="50%" style="display: block; margin: auto;" /> --- # Regresión lineal con error de tipo ARIMA ```r (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)$$` --- # Regresión lineal con error de tipo ARIMA <img src="presentacion_files/figure-html/unnamed-chunk-5-1.png" width="50%" style="display: block; margin: auto;" /> --- # Regresión lineal con error de tipo ARIMA .pull-left[ - Residuales del modelo de regresión con errores independientes. <img src="presentacion_files/figure-html/unnamed-chunk-6-1.png" width="90%" style="display: block; margin: auto;" /> ] .pull-right[ - Residuales del modelo dinámico con error de estructura ARMA. <img src="presentacion_files/figure-html/unnamed-chunk-7-1.png" width="90%" style="display: block; margin: auto;" /> ] --- # Regresión lineal con error de tipo ARIMA <img src="presentacion_files/figure-html/unnamed-chunk-8-1.png" width="40%" style="display: block; margin: auto;" /> ``` ## ## 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 ``` --- # Regresión lineal con error de tipo ARIMA ## Pronóstico: ```r fcast <- forecast(mod, xreg=rep(mean(uschange[,2]),8)) autoplot(fcast) + xlab("Year") + ylab("Percentage change") ``` <img src="presentacion_files/figure-html/unnamed-chunk-9-1.png" width="40%" style="display: block; margin: auto;" /> --- ## 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\)`. --- ## Serie estacionaria por tendencia y por diferencia. - 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. --- ## Serie estacionaria por tendencia y por diferencia. - 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: 1. tendencia determinística, 2. tendencia estocástica, y 3. 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. --- ## Serie estacionaria por tendencia y por diferencia. .pull-left[ <img src="presentacion_files/figure-html/unnamed-chunk-10-1.png" width="80%" style="display: block; margin: auto;" /> 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\)` ] .pull-right[ <img src="presentacion_files/figure-html/unnamed-chunk-11-1.png" width="80%" style="display: block; margin: auto;" /> 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. --- # Tendencia determinística y estocástica <img src="presentacion_files/figure-html/unnamed-chunk-12-1.png" width="50%" style="display: block; margin: auto;" /> --- # Tendencia determinística y estocástica ```r 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)$$` --- # Tendencia determinística y estocástica ```r (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)$$` --- # Tendencia determinística y estocástica <img src="presentacion_files/figure-html/unnamed-chunk-15-1.png" width="500px" height="300px" style="display: block; margin: auto;" /> --- # 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). --- # Análisis de intervención - 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. <img src="presentacion_files/figure-html/unnamed-chunk-16-1.png" width="400px" height="300px" style="display: block; margin: auto;" /> --- # Análisis de intervención **Pasos a seguir:** 1. Elaborar un modelo ARIMA antes de la intervención, i.e. durante las primeras 134 semanas. 2. Ampliar el modelo agregando variables indicadoras para registrar la intervención. 3. Re-estimar el modelo con la serie completa con las variables indicadoras. --- # Análisis de intervención .pull-left[ <img src="presentacion_files/figure-html/unnamed-chunk-17-1.png" width="90%" style="display: block; margin: auto;" /> ] .pull-right[ - Crest1: antes de la intervención <img src="presentacion_files/figure-html/unnamed-chunk-18-1.png" width="90%" /> ] --- # Análisis de intervención - Serie diferenciada (antes de la intervención) .pull-left[ <img src="presentacion_files/figure-html/unnamed-chunk-19-1.png" width="90%" style="display: block; margin: auto;" /> ] .pull-right[ <img src="presentacion_files/figure-html/unnamed-chunk-20-1.png" width="90%" /> ] --- # Análisis de intervención - Vamos a ajustar un ARIMA(0,1,1): `$$(1-B)Z_t=(1-\theta_1B)a_t$$` ```r 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 ``` --- # Análisis de intervención - 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. $$ $$ I_{2t}=\left\lbrace `\begin{aligned} 0, & & \text{si}~~ t<136 \\ 1, & & \text{si}~~ t \geq136, \end{aligned}` \right. $$ --- # Análisis de intervención - 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$$` --- # Análisis de intervención ```r 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) 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 ``` --- # Análisis de intervención <img src="presentacion_files/figure-html/unnamed-chunk-23-1.png" width="40%" style="display: block; margin: auto;" /> ``` ## ## 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 ```