Modelos Multivariados
Los datos se suelen recolectar en base a más de una variable. Por ejemplo, en economía, las tasas de interés diarias están disponibles para un amplio rango de monedas. Es por eso que es necesario realizar modelos donde no sólo se tome a una variable. En este tópico se introducirán modelos en los que se analiza la relación de diversas variables de series de tiempo. Se ahondará en los temas correspondientes al análisis de raíces unitarias en serie temporales, que servirá para evaluar la estacionariedad de la serie.
Regresión Espuria
Es una práctica común utilizar la regresión para explorar la relación entre dos o más variables, y usualmente se buscan predictores que tengan que ver directamente con la variable de respuesta y dan una explicación razonable para la variable dependiente. Para las series de tiempo se debe ser particularmente cuidadosos y sigilosos antes de introducir una relación causal, ya que esta aparente relación puede existir por el simple hecho de la tendencia que presentan las series y no por la relación en sí misma. A continuación, se presenta un gráfico de dispersión entre el producto en millones de soles del sector pesquero y del sector construcción y el coeficiente de correlación asociado a dichas variables.Después de subir las series de datos, se transforman en series de tiempo y luego se grafican.
pesca.ts = ts(data_R_ava_6$pesca, start = 1950, freq = 1)
const.ts = ts(data_R_ava_6$const, start = 1950, freq = 1)
plot(as.vector(aggregate(pesca.ts)), as.vector(aggregate(const.ts)))
El coeficiente de correlación sería de la siguiente manera:
## [1] 0.8267451La alta correlación de 0.83 y el gráfico de dispersión muestran una alta correlación de forma directa entre el producto del sector pesca y el del sector construcción; pero estos resultados no implican que exista una relación causal entre ambas variables. Se puede ajustar un modelo de regresión de una variable como función lineal de la otra, pero esta regresión va a ser llamada espuria por la relación causal inexistente entre las variables. El término regresión espuria es usado también cuando la tendencia estocástica subyacente en ambas series coincide. Las tendencias estocásticas son características de un proceso 𝐴𝑅𝐼𝑀𝐴 con raíz unitaria (𝐵 = 1 es una solución de la ecuación característica). Esto puede ser ilustrado simulando dos caminos aleatorios:
set.seed(10); x <- rnorm(100); y <- rnorm(100)
for(i in 2:100) {
x[i] <- x[i-1] + rnorm(1)
y[i] <- y[i-1] + rnorm(1) }
plot(x, y)
El coeficiente de correlación se estima de la siguiente manera:
## [1] 0.9039082Proceso de Raíz Unitaria
Se tiene el siguiente proceso AR(1):
\[ x_t = px_{t-1} + w_t, \text{ } -1 \leq p \leq 1 \]
Donde \(w_t\) es un ruido blanco.
Un proceso de raíz unitaria se va a dar cuando p = 1, lo cual llevará a que dicho proceso sea un camino aleatorio puro, y se tendría por lo tanto un proceso no estacionario para \(x_t\). Si a la expresión anterior se la resta \(x_{t-1}\) en ambos lados, se tendrá que:
\[ \vartriangle x_t = \delta x_{t-1} + w_t \]
Donde \(\delta = p - 1\). Así que en la práctica se estimará el modelo de la última expresión, para probar la hipótesis nula de que \(\delta = 0\); es decir, existe raíz unitaria y por lo tanto la serie \(x_t\) no es estacionaria. El estadístico para realizar dicha hipótesis sigue una distribución del estadístico \(tau (\tau)\).
Pruebas de Raíz Unitaria
Dentro de las pruebas que se tiene para probar la existencia de raíz unitaria se tiene la prueba de Dickey Fuller (DF), Dickey Fuller Aumentado (ADF), Phillips-Perron (PP) y Kwiatkowski-Phillips-Smichdt-Shin (KPSS).
Dickey Fuller
Es el modelo más simple para evaluar la presencia de raíz unitaria, este modelo viene dado por:
\[ \vartriangle x_t = \alpha + \delta x_{t-1} + w_t \]
El contraste será el siguiente:
\(H_0 \delta = 0 \rightarrow\) Existe raíz unitaria, \(x_t\) entonces la serie no es estacionaria.
\(H_0 \delta \neq 0 \rightarrow\) No existe raíz unitaria, \(x_t\) entonces la serie es estacionaria.
Si:
\(\tau - calculado\) en valor absoluto > \(\tau - crítico\) en valor absoluto: Se rechaza la \(H_o\).
\(\tau - calculado\) en valor absoluto < \(\tau - crítico\) en valor absoluto: Se acepta la \(H_o\).
Para realizar la prueba de Dickey-Fuller en R se hace uso de la función ur.df() de la librería urca.
El gráfico de la serie es el siguiente:
Se hará el análisis de raíz unitaria:
Los resultados podrán verso haciendo un summary:
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.31407 -0.86346 0.07963 0.66540 2.03542
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.009083 0.031521 -0.288 0.774
##
## Residual standard error: 0.9757 on 98 degrees of freedom
## Multiple R-squared: 0.0008466, Adjusted R-squared: -0.009349
## F-statistic: 0.08304 on 1 and 98 DF, p-value: 0.7738
##
##
## Value of test-statistic is: -0.2882
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.6 -1.95 -1.61Ya que el valor calculado de 0.288 es menor en términos absolutos a lo valores de tau, la hipótesis nula es aceptada, por lo tanto, existe raíz unitaria y la serie no es estacionaria. Otra forma es haciendo uso de la función adf.test() de la librería tseries.
## Warning: package 'tseries' was built under R version 4.0.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Augmented Dickey-Fuller Test
##
## data: x
## Dickey-Fuller = -1.9334, Lag order = 0, p-value = 0.6042
## alternative hypothesis: stationaryEl p-valor es 0.6042, lo que nos dice que la hipótesis nula es aceptada y por lo tanto la serie no es estacionaria.
Dickey Fuller Aumentado
En esta prueba se puede excluir la constante e incluir una tendencia lineal. La prueba ADF consiste en la estimación del siguiente modelo:
\[ \vartriangle x_t = \beta_0 + \beta_1 t + \delta x_{t-1} + \alpha_i \sum_{i=1}^m \vartriangle x_{t-i} + w_t \]
El contraste es similar al caso de la prueba Dickey-Fuller:
• 𝐻0: 𝛿 = 0 → Existe raíz unitaria, 𝑥𝑡 no es estacionaria.
• 𝐻1: 𝛿 ≠ 0 → No existe raíz unitaria, 𝑥𝑡 es estacionaria.
Si:
• 𝜏 -calculado en valor absoluto > 𝜏 -crítico en valor absoluto: Se rechaza 𝐻0.
• 𝜏 -calculado en valor absoluto < 𝜏 -crítico en valor absoluto: Se acepta 𝐻0.
Para hacer la estimación de dicho test en R se usan las mismas funciones que para el Dickey-Fuller simple, pero especificando en el número de rezagos igual a 1. Se puede observar que en ambos casos la hipótesis nula es aceptada y esto da razón de la existencia de raíz unitaria en la serie, por lo tanto, esta es estacionaria.
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.25453 -0.90101 0.02681 0.70255 1.97744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.01413 0.03282 -0.430 0.668
## z.diff.lag 0.06424 0.10583 0.607 0.545
##
## Residual standard error: 0.983 on 96 degrees of freedom
## Multiple R-squared: 0.00467, Adjusted R-squared: -0.01607
## F-statistic: 0.2252 on 2 and 96 DF, p-value: 0.7988
##
##
## Value of test-statistic is: -0.4304
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.6 -1.95 -1.61Phillips Pherron
Esta prueba estima una regresión haciendo una corrección sobre lamatriz de varianzas y covarianzas de los residuos. La corrección esmediante un método no paramétrico. En esta se estima la siguiente regresión:
\[ \vartriangle x_t = \vartriangle \beta + px_{t-1} + w_t \]
A diferencia de la prueba ADF, no existen términos de diferencias retardados. La hipótesis nula 𝐻0 del test de Phillips-Perron es la trayectoria de raíz unitaria con tendencia y la alternativa la estacionariedad con tendencia, si el valor 𝑡 asociado al coeficiente de 𝑥𝑡−1 es mayor en valor absoluto al estadístico de MacKinnon, se rechaza la hipótesis nula de existencia de raíz unitaria. test de Phillips-Perron en R, estas son PP.test, pp.test (librería tseries) y ur.pp (librería urca).
- PP.test:
##
## Phillips-Perron Unit Root Test
##
## data: x
## Dickey-Fuller = -2.1739, Truncation lag parameter = 3, p-value = 0.5046- pp.test:
##
## Phillips-Perron Unit Root Test
##
## data: x
## Dickey-Fuller Z(alpha) = -10.323, Truncation lag parameter = 3, p-value
## = 0.5169
## alternative hypothesis: stationary- ur.pp:
##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept
##
##
## Call:
## lm(formula = y ~ y.l1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2061 -0.7372 0.1559 0.7757 2.1076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10167 0.10128 -1.004 0.318
## y.l1 0.98274 0.03256 30.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9756 on 97 degrees of freedom
## Multiple R-squared: 0.9038, Adjusted R-squared: 0.9028
## F-statistic: 911.2 on 1 and 97 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic, type: Z-alpha is: -2.6159
##
## aux. Z statistics
## Z-tau-mu -0.9923En los tres se observa que la hipótesis nula se acepta y se entenderá que la serie no es estacionaria.
Kwiatkowski-Phillips-Smichdt-Shin (KPSS)
En esta prueba se propone contrastar como hipótesis nula la estacionariedad de las tendencias, esta es la principal diferencia con las pruebas de raíz unitaria.
El contraste es el siguiente:
• 𝐻0: La serie es estacionaria. • 𝐻1: La serie no es estacionaria.
Si:
• Valor calculado > Valor crítico: Se rechaza la 𝐻0. • Valor calculado < Valor crítico: Se acepta la 𝐻0.
Para realizar en R el contraste KPSS, se debe hacer uso de las funciones ur.kpss de la librería urca o de la función kpss.test.
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 4 lags.
##
## Value of test-statistic is: 1.2788
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739Se puede observar para ambos casos que la hipótesis nula es rechazada, por lo tanto, la serie no es estacionaria, tal y como se había hallado con las pruebas anteriores.
Cointegracíon
Existen modelos donde la relación de series no es espuria, es decir, ambas series tendrán una correlación alta y esta será producto de la relación directa entre ambas. Muchas series están altamente correlacionadas en el tiempo. Por ejemplo, se puede observar que las tasas de interés de la libra en el Reino Unido y el Euro muestran una alta correlación. Esto puede ser explicado por la similitud de estas economías con la de Estados Unidos. En el anterior ejemplo respecto a la relación del producto entre el sector pesca y el sector construcción en Perú, que puede ser explicado por el aumento de un mismo factor para ambas, pero que estas dos no tienen relación de forma directa o causal. Se puede observar también que dos series que son independientes y contienen raíces unitarias, pueden mostrar una aparente relación lineal, la relación entre ambas series recibe el nombre de espuria. Sin embargo, se puede observar que dos series pueden contener raíces unitarias y estar relacionadas, cuando sucede esto se dice que estas series están cointegradas. La definición estricta de cointegración es la siguiente:
Dos series estacionarias \({x_t}\) y \({y_t}\) son cointegradas si alguna combinación lineal \(ax_t+by_t\) con constante a y b, es una serie estacionaria. Como ejemplo se puede considerar un camino aleatorio 𝜇𝑡, dado por 𝜇𝑡 = 𝜇𝑡−1 + 𝑤𝑡, donde 𝑤𝑡 es un ruido blanco con media cero, y dos series 𝑥𝑡 & 𝑦𝑡, dadas por 𝑥𝑡 = 𝜇𝑡 + 𝑤𝑥,𝑡 y 𝑦𝑡 = 𝜇𝑡 + 𝑤𝑦,𝑡, donde 𝑤𝑥,𝑡 y 𝑤𝑦,𝑡 son series de ruido blanco independiente con media cero. Ambas series son no estacionarias, pero su diferencia 𝑥𝑡 − 𝑦𝑡 es estacionaria, ya que es una combinación lineal de los términos de ruido blanco independiente. Así, la combinación lineal de 𝑥𝑡 y 𝑦𝑡 , con 𝑎 = 1 y 𝑏 = −1 , produce una serie estacionaria. Por lo tanto 𝑥𝑡 y 𝑦𝑡 están cointegradas y comparten la tendencia estocástica subyacente 𝜇𝑡. En R, dos series pueden ser probadas para cointegración haciendo uso del test de Phillips-Oularis, implementado en la función po.test que está dentro de la librería tseries. La función requiere que las series estén dadas en forma matricial, este test produce los resultados para un test donde la hipótesis nula nos dice que las dos series no están cointegradas. Como ejemplo, se simula dos series cointegradas que tienen una tendencia estocástica y se usa el comando po.test para probar que estas series son cointegradas:
x <- y <- mu <- rep(0, 1000)
for (i in 2:1000) mu[i] <- mu[i - 1] + rnorm(1)
x <- mu + rnorm(1000)
y <- mu + rnorm(1000)Se realiza una prueba de raíz unitaria para ambas series:
## [1] 0.574204## [1] 0.6089259Ambas series presentan raíces unitarias, por lo que no son estacionarias, la sintaxis para realizar el test de cointegración es la siguiente:
## Warning in po.test(cbind(x, y)): p-value smaller than printed p-value##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(x, y)
## Phillips-Ouliaris demeaned = -1034.1, Truncation lag parameter = 9,
## p-value = 0.01Los resultados del test estimado nos dan evidencia de que las series están cointegradas ya que la hipótesis nula es rechazada al 1% de significancia.
Modelo de Vectores Autoregresivos
Dos series de tiempo, 𝑥𝑡 y 𝑦𝑡 , siguen un proceso de vectores autorregresivos de orden 1 (denotado como 𝑉𝐴𝑅(1)) si:
\[ x_t = \theta_{11} x_{t-1} + \theta_{12} y_{t-1} + w_{x,t} y_t = \theta_{21} x_{t-1} + \theta_{22} t_{t-1} + w_{y,t} \]
Donde 𝑤𝑥,𝑡 y 𝑤𝑦,𝑡 son ruidos blancos bivariados y 𝜃𝑖𝑗 son los parámetros del modelo. Si las secuencias de ruido blanco están definidas con media 0 y el proceso es estacionario, ambas series, 𝑥𝑡 y 𝑦𝑡 tienen media 0. El sistema definido anteriormente puede ser reescrito en notación matricial como:
\[ Z_t = \theta Z_{t-1} + w_t \]
Donde:
La última ecuación es una expresión en vectores de un proceso 𝐴𝑅(1), es decir, el proceso es un vector autorregresivo. Haciendo uso del operador de rezago, la ecuación puede ser reescrita como:
\[ (I-\theta B) Z_t = \theta (B) Z_t = w_t \]
Donde 𝜃 es una matriz polinomial de orden 1 e 𝐼 es una matriz identidad 2×2. Un proceso 𝑉𝐴𝑅(1) puede ser extendido a un 𝑉𝐴𝑅(𝑝) permitiendo a 𝜃 ser una matriz polinomial de orden 𝑝. Un modelo 𝑉𝐴𝑅(𝑝) para 𝑚 series de tiempo está definido en la ecuación anterior, donde 𝐼 es la matriz identidad 𝑚×𝑚, 𝜃 es un polinomio de matrices 𝑚×𝑚 de parámetros, 𝑍𝑡 es una matriz 𝑚×1 de variables de series de tiempo, y 𝑤𝑡 es el ruido blanco multivariado. Para un modelo 𝑉𝐴𝑅, la ecuación característica está dada por el determinante de una matriz. Análogamente a los modelos 𝐴𝑅, un 𝑉𝐴𝑅(𝑝) es estacionario si las raíces del determinante |𝜃(𝑥)| exceden todas a la unidad en valor absoluto. Para el modelo 𝑉𝐴𝑅(1), el determinante está dado por:
La función de R polyroot y Mod pueden ser usadas para probar si un modelo 𝑉𝐴𝑅 es estacionario, donde la función polyroot toma un vector de coeficientes polinómicos como entrada del parámetro. Por ejemplo, considerando a un modelo 𝑉𝐴𝑅(1) con una matriz de parámetros \(\theta = 0.4, 0.3, 0.2, 0.1\). Entonces la ecuación característica estará dada por:
El valor absoluto de las raíces de la ecuación está dado por:
## [1] 1.861407 26.861407De esto se puede deducir que el modelo 𝑉𝐴𝑅(1) es estacionario ya que ambas raíces exceden a la unidad en valor absoluto. Los parámetros de un modelo 𝑉𝐴𝑅(𝑝) pueden ser estimados usando la función ar en R, que selecciona el mejor ajuste del orden 𝑝 basado en el 𝐴𝐼𝐶 más pequeño. Usando el proceso bivariado de proceso de ruido blanco de la ecuación (6) y los parámetros ya dados del proceso 𝑉𝐴𝑅(1), un proceso 𝑉𝐴𝑅(1) es simulado y los parámetros son estimados usando la función ar:
set.seed(1)
x <- y <- rep(0, 1000)
wx = rnorm(1000); wy = rnorm(1000)
x[1] <- wx[1]
y[1] <- wy[1]
for (i in 2:1000) {
x[i] <- 0.4 * x[i - 1] + 0.3 * y[i - 1] + wx[i]
y[i] <- 0.2 * x[i - 1] + 0.1 * y[i - 1] + wy[i]
}
xy.ar <- ar(cbind(x, y))
xy.ar$ar[, , ]## x y
## x 0.3474962 0.2827365
## y 0.2314242 0.1146880Como era de esperarse, los parámetros estimados son cercanos a los valores subyacentes del modelo.
Modelo VAR para las series económicas de EE.UU
Una serie de datos trimestrales (1954-1987), están disponibles en la librería tseries. El modelo 𝑉𝐴𝑅 de mejor ajuste estima los parámetros de un modelo con el producto nacional bruto (GNP) y la cantidad total de dinero (M1) en el siguiente ejemplo. Mínimos cuadrados ordinarios son usados para ajustar el modelo.
Ejemplos
Ejemplo 1: Cointegración
La base de datos siguiente se muestra información sobre la evolución de la producción de espárragos en miles de toneladas y el índice del producto bruto interno de forma mensual, desde enero de 2007 hasta noviembre de 2016. En el siguiente ejemplo se estimará un modelo de regresión para ambas variables y se probará que existe cointegración en el mismo haciendo uso de R como herramienta principal.
Se declaran las series de tiempo:
pbi.ts = ts(pbi_esparr$pbi, start=2007, freq = 12)
esparr.ts = ts(pbi_esparr$esparr, start=2007, freq = 12)Podemos hacer un gráfico de ambas series con la función ts.plot():
Se puede observar una ligera tendencia al alza en las dos series, se procede a realizar el test de cointegración de Phillips-Ouliaris:
## Warning in po.test(cbind(pbi.ts, esparr.ts)): p-value greater than printed p-
## value##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(pbi.ts, esparr.ts)
## Phillips-Ouliaris demeaned = -11.776, Truncation lag parameter = 1,
## p-value = 0.15La hipótesis nula es aceptada, y se llega a la conclusión de que las series en cuestión no están cointegradas.
Ejemplo 2: VAR
La base de datos siguiente muestra información sobre dos variables la demanda interna (DI) y la cantidad real de dinero (MP). En este ejemplo se estimará un VAR(2) y se realizarán predicciones para ambas series a partir de dicho modelo estimado.
Podemos hacer un gráfico de ambas series con la función ts.plot():
Se puede observar una tendencia al alza en ambas series. Para realizar la estimación del modelo VAR(2) cargamos en primer lugar la librería vars, luego de eso procedemos con la estimación del modelo:
## Warning: package 'vars' was built under R version 4.0.5## Loading required package: MASS## Warning: package 'MASS' was built under R version 4.0.4## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 4.0.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.0.2## Loading required package: lmtestPara ver los coeficientes seguimos la siguiente sintaxis:
## $di
## Estimate Std. Error t value Pr(>|t|)
## di.l1 0.2424892 8.997126e-02 2.695185 7.929673e-03
## mp.l1 1.1289965 3.656866e-01 3.087334 2.451649e-03
## di.l2 0.6056329 8.970477e-02 6.751402 3.963445e-10
## mp.l2 -0.6179203 3.718730e-01 -1.661643 9.890480e-02
## const 3419.8965270 2.183512e+03 1.566236 1.196344e-01
##
## $mp
## Estimate Std. Error t value Pr(>|t|)
## di.l1 1.872871e-02 0.02594140 0.72196217 4.715664e-01
## mp.l1 9.427476e-01 0.10543836 8.94122056 2.575799e-15
## di.l2 1.036235e-03 0.02586456 0.04006390 9.681014e-01
## mp.l2 6.775089e-03 0.10722210 0.06318743 9.497107e-01
## const -3.692753e+02 629.57184889 -0.58654985 5.584856e-01Luego de esto se puede generar predicciones para el modelo estimado, haciendo uso de la función predict():
Se están generando pronósticos para los ocho periodos posteriores, luego de esto se genera una serie de tiempo con los pronósticos estimados:
di.pred <- ts(var.pred$fcst$di[, 1], st = c(2015,3), fr = 4)
mp.pred <- ts(var.pred$fcst$mp[, 1], st = c(2015,3), fr = 4)Luego de haber generado la serie de pronósticos veremos que se puede observar en un gráfico la serie:
