UNIVERSIDAD DE EL SALVADOR

FACULTAD DE CIENCIAS ECONÓMICAS

ESCUELA DE ECONOMÍA

CICLO II-2021

Asignatura:

Métodos para el Análisis Ecocómico.

Facilitador:

Carlos Ademir Pérez Alas.

Contenido:

Laboratorio 3.

Grupo Teórico:

03.

Integrantes:

Hernández Cartagena, Karla Mireya. HC18035.

Mejía Ramos, Carlos Arnoldo. MR18006.

Meléndez Morales, Kennya Elizabeth. MM16049.

Novoa Rubio, Fernando Enrique. NR18001.

Ciudad Universitaria, viernes 3 de diciembre de 2021.

Laboratorio 3

Prueba de raíz unitaria de Dickey & Fuller

Propósito

La Prueba Dickey Fuller es una prueba que permite determinar la presencia de raíz unitaria en series de tiempo, a partir de esto se demuestra la existencia de estacionariedad o no estacionariedad. Esta prueba es una de las más utilizadas para este fin, siendo importante porque para algunos modelos como los modelos ARIMA se requiere que la serie de tiempo sea estacionaria. (Rico, n.d.)

Hipótesis

Si consideramos un proceso estocástico de la forma:

\[Y_{t} = \phi Y_{t-1} + u_{t}\ , \ -1 \leq \phi \leq 1\]

donde: \(Y_{t}\) es la serie de tiempo y \(u_{t}\) es un término de error de ruido blanco.

\(H_{0}:\) \(\beta = 0\) (es equivalente a \(\phi\) = 1)

“Es decir, existe raiz unitaria; por lo tanto, no hay estacionariedad”

\(H_{1}:\) \(\beta < 0\) (es equivalente a \(\phi\) < 1)

“Es decir, no existe raiz unitaria; por lo tanto, hay estacionariedad”

Sintaxis

Librerías a utilizar:

library(urca)

Prueba de Dickey-Fuller

ur.df (y, type = c ("none", "drift", "trend"), lags = 1, selectlags = c ("Fixed", "AIC", "BIC"))

Argumentos:

y: Es el vector que se probará para la raiz unitaria (residuos).
type: Es el tipo de prueba (“none” ni una intersección ni una tendencia, “drift” una intersección, “trend” una intersección y una tendencia).
selectlags: Selección de los retardos. Esta se puede lograr de acuerdo con los criterios de Akaike “AIC” o Bayes “BIC.” El número máximo de rezagos considerados lo establece ´lags´. El valor predeterminado es utilizar una “Fixed”longitud de retraso establecida por ´lags´.

Prueba de Dickey Fuller Aumentada

adf.test (x, nlag = NULL, output = TRUE)

Argumentos:

x = un vector numérico o una serie de tiempo univariante.
nlag = el orden de retraso con el valor predeterminado para calcular la estadística de prueba.
output = un valor lógico que indica imprimir los resultados de la prueba en la consola R. El valor predeterminado es TRUE.

Estadístico de Prueba

Estadístico \(\tau(Tau)\) sobre la variable de prueba rezagada, se utiliza \(\tau\) debido a que las pruebas \(t\) usuales no son adecuadas. La salida que R nos da contiene un valor t que es el que contrastaremos con el valor crítico,

Criterio de decisión

\(No\ se\ Rechaza\ H_{0}: \beta=0 \ ó \ \ P-Value=1\)

La serie de tiempo no es estancionaria es decir tiene tendencia estocástica.

\(Se\ rechaza\ H_{0}: \beta<0 \ ó \ \ P-Value<1\)

La serie de tiempo es estancionaria es decir tiene una tendencia determinista.

Interpretación del Estadístico de Prueba

Si el estadístico de prueba es mayor que el valor critico se tendrá la evidencia suficiente en favor de no rechazar la hipotesis nula es decir que la serie no es estacionaria.

Si el estadístico de prueba es menor que el valor critico se tendrá la evidencia suficiente en favor de rechazar la hipotesis nula es decir que la serie es estacionaria.

Ejemplo

El ejemplo consiste en saber si las series de tiempo DPI (Ingreso personal disponible) y PCE (Gasto de consumo personal) de la base de datos coint son o no son estacionarias. (Cuellar 2020)

Para cargar la data es necesaria la libreria readxl (Wickham and Bryan 2019), para el formato de tabla “kableExtra” (Zhu 2021) “car” carga la función residualPlot (Fox and Weisberg 2019)

library(tseries) # carga la función adf.test
library(urca)
library(car) # carga la función residualPlot
library(readxl)
library(kableExtra)
coint <- read_excel("coint.xls")
kable(
  head(coint, 10),
  caption = "**Tabla 1:** Muestra de 10 Observaciones de la serie de tiempo coint",
  align = 'c',
  row.names = TRUE,
  digits = 2
) %>%
  kable_classic(html_font = "Times New Roman",
                font_size = 14) %>%
  row_spec(0, bold = T) %>%
  footnote(general_title = "**Fuente**:",
           general = "Elaboración propia con base en datos proporcionados por Cuellar 2020")

**Tabla 1:** Muestra de 10 Observaciones de la serie de tiempo coint
	DPI	PCE
1	1096.0	1017.2
2	1072.8	1034.0
3	1102.8	1037.5
4	1089.7	1037.7
5	1107.3	1042.6
6	1145.3	1054.3
7	1168.4	1056.1
8	1171.9	1064.8
9	1147.6	1066.1
10	1151.4	1082.6
Fuente:
Elaboración propia con base en datos proporcionados por Cuellar 2020

# Generar Logaritmos  
lnPCE = log(coint$PCE)
lnDPI = log(coint$DPI)

# Generar Variables de Series de Tiempo
DPI.ts <- ts(lnDPI, start = c(1947, 1), frequency = 4)
PCE.ts <- ts(lnPCE, start = c(1947, 1), frequency = 4)

# Hacer una tabla
datos1 <- cbind(DPI.ts, PCE.ts)
plot(datos1, main = "Tendencia")

# Generar Modelo
modelo1 <- lm(PCE.ts ~ DPI.ts)
summary(modelo1)

## 
## Call:
## lm(formula = PCE.ts ~ DPI.ts)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.062816 -0.019375 -0.001032  0.017186  0.078186 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -0.194233   0.023593  -8.233   0.0000000000000114 ***
## DPI.ts       1.011351   0.002902 348.542 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02764 on 242 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.998 
## F-statistic: 1.215e+05 on 1 and 242 DF,  p-value: < 0.00000000000000022

residuales <- modelo1$residuals
residualPlot(modelo1)

# Prueba de Dickey-Fuller
y <- ur.df(residuales)
summary(y)

## 
## ############################################### 
## # 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 
## -0.050339 -0.004034 -0.000311  0.005048  0.042822 
## 
## Coefficients:
##            Estimate Std. Error t value    Pr(>|t|)    
## z.lag.1    -0.05978    0.02364  -2.529      0.0121 *  
## z.diff.lag -0.31809    0.05924  -5.369 0.000000187 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009841 on 240 degrees of freedom
## Multiple R-squared:  0.1503, Adjusted R-squared:  0.1433 
## F-statistic: 21.23 on 2 and 240 DF,  p-value: 0.000000003228
## 
## 
## Value of test-statistic is: -2.5289 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

y@teststat

##                tau1
## statistic -2.528947

y@cval

##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Para utilizar la funcion ur.df se usa la libreria “urca” (Pfaff 2008)

#Prueba con intercepto
y2 <- ur.df(residuales, type = "drift", selectlags = "AIC")
summary(y2)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.050207 -0.003904 -0.000181  0.005176  0.042952 
## 
## Coefficients:
##               Estimate Std. Error t value    Pr(>|t|)    
## (Intercept) -0.0001294  0.0006339  -0.204      0.8385    
## z.lag.1     -0.0598500  0.0236882  -2.527      0.0122 *  
## z.diff.lag  -0.3180310  0.0593632  -5.357 0.000000198 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00986 on 239 degrees of freedom
## Multiple R-squared:  0.1504, Adjusted R-squared:  0.1433 
## F-statistic: 21.15 on 2 and 239 DF,  p-value: 0.000000003478
## 
## 
## Value of test-statistic is: -2.5266 3.2058 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

Esta prueba considera la tendencia e intecepto, mostrando un valor estadístico de \(-2.53\) es menor en valor absoluto a valor critico al \(10\%\) de \(-2.57\), comprobando la existencia de raíz unitaria por lo que se concluye que el PCI y el DPI no tienen una relación de cointegración.

#Prueba sin tendencia y sin intercepto
y3 <- ur.df(residuales, type = "none", selectlags = "AIC")
summary(y3)

## 
## ############################################### 
## # 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 
## -0.050339 -0.004034 -0.000311  0.005048  0.042822 
## 
## Coefficients:
##            Estimate Std. Error t value    Pr(>|t|)    
## z.lag.1    -0.05978    0.02364  -2.529      0.0121 *  
## z.diff.lag -0.31809    0.05924  -5.369 0.000000187 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009841 on 240 degrees of freedom
## Multiple R-squared:  0.1503, Adjusted R-squared:  0.1433 
## F-statistic: 21.23 on 2 and 240 DF,  p-value: 0.000000003228
## 
## 
## Value of test-statistic is: -2.5289 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

# lc.df <- ur.df(
#   residuales,
#   type = 'trend',
#   lags = 4,
#   selectlags = c("AIC")
# )
# summary(lc.df)

Esta prueba no considera el intecepto, mostrando un valor estadístico de \(-2.53\) es mayor en valor absoluto a valor critico al \(10\%\) de \(-1.62\), por lo que se concluye que el PCI y el DPI tienen una relación de cointegración.

#Dickey Fuller Aumentada
adf.test(residuales)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  residuales
## Dickey-Fuller = -1.2524, Lag order = 6, p-value = 0.8904
## alternative hypothesis: stationary

El chunk anterior empleo la libreria “tseries” (Trapletti and Hornik 2021)

Cointegración en el enfoque de Soren Johansen.

Propósito

La prueba de Soren Johansen se utiliza para verificar la existencia o no de cointegraccion entre dos o mas series temporales. (Mendoza González and Quintana Romero, n.d.)

Hipótesis

\(H_{0}: \ r=0 \ No \ existe \ cointegración\)

\(H_{1}: \ r=m \ Existe \ cointegración\)

Sintaxis

Librerías a utilizar:

library(urca)

Dentro de esta libreria se encuentra la función “ca.jo” que brinda el cálculo de la cointegración por el enfoque de johansen.

Prueba de Soren Johansen

johansen <- ca.jo(x, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2,
spec=c("longrun", "transitory"), season = NULL, dumvar = NULL)

summary(johansen)

Argumentos:

X: Matriz de datos a investigar para la cointegración.
type: La prueba a realizar, ya sea eigen o trace.
ecdet: Carácter, none para no interceptar en la cointegración, const por término constante en cointegración y trend para la variable de tendencia en cointegración
k: El orden de retraso de la serie (niveles) en el VAR.
spec:Determina la especificación del VECM, consulte los detalles a continuación.
season: Si se deben incluir variables ficticias estacionales, la frecuencia de los datos debe establecerse en consecuencia, es decir 4 para datos trimestrales.
dumvar: Si se deben incluir variables ficticias, se puede proporcionar una matriz con una dimensión de fila igual a x.

El summary es para la visualización de la prueba

Estadístico de Prueba

Para la salida de ambas pruebas, los estadisticos de estas se ubicaran en la columna test del resumen de estadisticas de pruebas y valores criticos.

Criterio de decisión

\(No\ se\ Rechaza\ H_{0}: \ r=0\)

En la serie de tiempo no existe cointegración.

\(Se\ Rechaza\ H_{0}: \ r=m\)

En la serie de tiempo existe cointegración.

El criterio de decisión para rechazar la hipótesis nula consiste el comparar el valor del estadistico de prueba y un valor critico para un determinado nivel de significancia, si el estadistico es mayor se rechaza la hipótesis nula, caso contrario no se rechaza la hipótesis nula.

Interpretación del Estadístico de Prueba

Si el estadístico de prueba es mayor que el valor critico se tendrán las suficientes evidencias para no aceptar la hipotesis nula es decir que en la serie de tiempo existe cointegración.

Ejemplo

Para realizar esta prueba se utilizara la misma data que en la prueba anterior lo primero que debemos hacer es cargar la data y llamar la libreia urca y car, luego nos aseguramos de sacar los logaritmos de las series temporales, una vez hecho esto procedemos a combinar las series en un mismo objeto y ya se puede aplicar la prueba.

A continuación la Prueba Johansen de cointegración de la traza

johansen <-
  ca.jo(
  coint,
  type = "trace",
  spec = "transitory",
  ecdet = "none",
  K = 2
  )
  summary(johansen)

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend 
## 
## Eigenvalues (lambda):
## [1] 0.183987546 0.007163676
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  1.74  6.50  8.18 11.65
## r = 0  | 50.94 15.66 17.95 23.52
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##           DPI.l1    PCE.l1
## DPI.l1 1.0000000  1.000000
## PCE.l1 0.7990846 -1.056443
## 
## Weights W:
## (This is the loading matrix)
## 
##            DPI.l1        PCE.l1
## DPI.d 0.002339272 -0.0243695001
## PCE.d 0.003349821  0.0008561397

Los resultados muestran que para la primera hipótesis \(r = 0\) y alternativa \(r = 1\), se rechaza la hipótesis nula de cointegración: el estadístico es mayor al valor critico del \(10\%\), \(50.94 > 15.66\).

Esta prueba de la traza se puede hacer con constante y con tendencia solo cambiando el valor de ecdet, \(ecdet = "const"\) para que sea constante y \(ecdet = "trend"\) para que sea con tendencia.

Prueba con constante

summary(ca.jo(
  coint,
  type = "trace",
  ecdet = "const",
  spec = c("longrun"),
  K = 4
  ))

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 0.118100383836122202474428 0.018683808020589515869192
## [3] 0.000000000000000001335927
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  4.53  7.52  9.24 12.97
## r = 0  | 34.69 17.85 19.96 24.60
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              DPI.l4       PCE.l4   constant
## DPI.l4     1.000000    1.0000000   1.000000
## PCE.l4    -0.214827   -0.9722601  -1.101688
## constant 716.842709 -499.9065251 -10.942805
## 
## Weights W:
## (This is the loading matrix)
## 
##            DPI.l4       PCE.l4                 constant
## DPI.d 0.005018441 -0.019139461  0.000000000000007891529
## PCE.d 0.004268926  0.004884931 -0.000000000000010101061

Prueba con tendencia

summary(ca.jo(
  coint,
  type = "trace",
  ecdet = "trend",
  spec = c("longrun"),
  K = 4
))

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1]  0.07018248829245908593233  0.04994817644718644450430
## [3] -0.00000000000000002775558
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 | 12.30 10.49 12.25 16.26
## r = 0  | 29.76 22.76 25.32 30.45
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              DPI.l4     PCE.l4 trend.l4
## DPI.l4    1.0000000  1.0000000  1.00000
## PCE.l4   -0.8691622 -0.8234953 -1.61412
## trend.l4 -8.9314092 -6.6891103 15.28862
## 
## Weights W:
## (This is the loading matrix)
## 
##            DPI.l4     PCE.l4                  trend.l4
## DPI.d -0.03848898 -0.1216981 -0.0000000000000294827677
## PCE.d -0.03374167  0.0312370  0.0000000000000005650615

Los resultados muestran que para la primera hipótesis \(r = 0\) y alternativa \(r = 1\), no se acepta la hipótesis nula de cointegración: el estadístico es mayor al valor critico del \(10%\), \(29.76 > 22.76\).

Causalidad en el sentido de Granger

Propósito

La prueba de causalidad de Granger se usa para examinar si una serie de tiempo puede usarse para pronosticar otra. (Lesmeister 2013) Con esta prueba pretendemos saber si la serie X causa a la serie Y o si la serie Y causa a la serie X (causalidad unidireccional), si hay una causalidad bilateral o si ninguna serie de tiempo causa a la otra.

Hipótesis

Como tenemos dos variables, la prueba implica la estimación de los siguientes modelos que predicen Y:

Uno con solo valores pasados de Y al que le denominamos \(\Omega\), donde k es el número de retardos de la serie temporal, este es denominado como modelo restringido.

\[ \Omega = Y_{t}= \beta_{0} + \beta_{1}Y_{t-1} + ... + \beta_{k}Y_{t-k} + e \]

Y otro con valores pasados de Y y X al que le denominamos \(\pi\), donde k es el número de retardos de la serie temporal, este es denominado como modelo no restringido.

\[ \pi = Y_{t}= \beta_{0} + \beta_{1}Y_{t-1} + ... + \beta_{k}Y_{t-k} + \alpha_{0} + \alpha_{1}X_{t-1} + ... + \alpha_{k}X_{t-k} + e \]

También resulta interesante realizar el modelo no restringido para la causalidad unidireccional de Y hacia X para verificar si existe una causalidad bilateral.

En este sentido las hipótesis son las siguientes:

Hipótesis nula:

\[ H_{0}: \alpha = 0, i = 1,2,...,n \]

Es decir, los términos rezagados de X no pertenecen a la regresión.

“Hay evidencia de que la serie de tiempo X (en el sentido de Granger) no causa la serie de tiempo Y.”

Hipótesis alterna:

\[H_{1}: \alpha_{i}\neq 0\ para\ al\ menos\ 1\ i\ del\ elemento\ [1, k]\]

Es decir, los términos rezagados de X pertenecen a la regresión.

“Hay evidencia de que la serie de tiempo X (en el sentido de Granger) causa la serie de tiempo Y.”

Esencialmente, estamos tratando de determinar con las hipótesis si X proporciona más información sobre los valores futuros de Y que únicamente los valores pasados de Y, por eso si los términos rezagados de X no pertenecen a la regresión, quiere decir que los movimientos en la serie de tiempo X no causa efecto en los valores de la serie de tiempo Y. (Lesmeister 2013)

Sintaxis

Para realizar la prueba de causalidad de Granger en R usamos la funcion grangertest de del paquete lmtest, su sintaxis es la siguiente:

grangertest(X, Y, order = 1)

Donde:

X: Es la primera serie temporal.

Y: Es la segunda serie temporal.

Order: Es el numero de retardos a utilizar en la primera serie de tiempo. El valor predeterminado es 1.(Lesmeister 2013)

Estadístico de Prueba

Para contrastar la hipótesis nula, aplicamos la prueba F dada por:

\[ F = \frac{(SCR_{R}-SCR_{NR})/m}{SCR_{NR}/(n-k)} \] Sigue la distribución F con m y (n - k) grados de libertad. En el presente caso, m es igual al número de términos rezagados de X, y k es el número de parámetros estimados en la regresión no restringida. (Gujarati and Porter 2010)

Ejemplo de la salida que genera R

library(lmtest)
grangertest(DAX ~ SMI, order = 3, data = tsData)
Granger causality test

Model 1: DAX ~ Lags(DAX, 1:3) + Lags(SMI, 1:3)
Model 2: DAX ~ Lags(DAX, 1:3)
  Res.Df Df      F    Pr(>F)    
1   1850                        
2   1853 -3 8.4968 1.322e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Al aplicar la prueba de causalidad de Granger podemos observar que el programa nos arroja los siguientes elementos:

Model 1: es el modelo no restringido que incluye los términos causales de Granger.
Model 2: es el modelo restringido donde se omiten los términos causales de Granger.

Abajo de esto en el numeral 3 nos muestra el resultado del valor del estadístico de prueba F y el valor del p-value con su respectivo nivel de significancia. Estos dos elementos son los que tomaremos en cuenta para el criterio de decisión sobre las hipótesis.

Criterio de decisión

Rechazar \(H_{0}\) si \(F \geq V.C\)
Rechazar \(H_{0}\) si \(P_{value} \leq \alpha\)

Se rechaza la hipótesis nula en caso de que el p-value sea menor que el nivel de significancia o el estadístico F sea mayor que el valor crítico. Como se observó en el ítem anterior, la salida de R nos proporciona estos dos elementos que utilizamos para contrastar las hipótesis de acuerdo al nivel de significancia elegido y el resultado del estadístico F crítico obtenido de la tabla de distribución F.

Interpretación

En caso de rechazar la hipótesis nula:

Existe evidencia suficiente a favor de rechazar la hipótesis nula, con lo que podemos concluir que los términos rezagados de X si pertenecen a la regresión y hay evidencia de que la serie de tiempo X (en el sentido de Granger) causa la serie de tiempo Y.

En caso de no rechazar la hipótesis nula:

Existe evidencia suficiente a favor de no rechazar la hipótesis nula, con lo que podemos concluir que los términos rezagados de X no pertenecen a la regresión y hay evidencia de que la serie de tiempo X (en el sentido de Granger) no causa la serie de tiempo Y.

Ejemplo de aplicación de la prueba

Para mostrar como se aplica la prueba de causalidad de Granger hemos tomado un ejemplo de (Lesmeister 2013) donde se utiliza la prueba para dar respuesta a laantigua pregunta ¿Qué fue primero, el huevo o la gallina?

Para dar respuesta ante este cuestionamiento se utilizan dos series de tiempo: la producción de huevos de Estados Unidos y la población estimada de pollos en Estados Unidos, ambas comprendidas desde 1930 hasta 1983. Estos datos fueron presentados por Walter Thurman y Mark Fisher en el American Journal of Agricultural Economics en mayo de 1988 con el título “Pollos, huevos y causalidad, o ¿cuál fue primero?” también vale la pena mencionar que Roger Koenker puso el set de datos ChickEgg a disposición de R aunque estos datos son ligeramente diferentes a los analizados por Thurman y Fisher. (Documentation, n.d.)

Carga de datos

library(dplyr)
library(kableExtra)
library(readxl)
ChickEgg <- read_excel("ChickEgg.xlsx")
ChickEgg %>% head() %>% kable(caption = "**Tabla:** Muestra de los primeros resultados de la población de pollos y producción de huevos en USA 1930 - 1983.",
align = "c",
digits = 6) %>% kable_classic(html_font = "Times New Roman",
font_size = 14) %>%
row_spec(0, bold = T) %>%
footnote(general_title = "**Fuente**:",
general = "Elaboración propia con base en datos de American Journal of Agricultural Economics.")

**Tabla:** Muestra de los primeros resultados de la población de pollos y producción de huevos en USA 1930 - 1983.
Year	chicken	egg
1930	468491	3581
1931	449743	3532
1932	436815	3327
1933	444523	3255
1934	433937	3156
1935	389958	3081
Fuente:
Elaboración propia con base en datos de American Journal of Agricultural Economics.

Para el uso del operador %>% es necesaria la libreria dplyr (Wickham et al. 2021) La serie de tiempo chicken muestra la población del 1 de diciembre de todos los pollos de Estados Unidos excluyendo los pollos de engorde comerciales.

La serie de tiempo egg muestra la producción de huevos de Estados Unidos en millones de docenas.

# Trazar series de tiempo
attach(ChickEgg)
par(mfrow = c(2, 1))
plot.ts(chicken)
plot.ts(egg)

Las gráficas proporcionan poca información aparte de que los datos probablemente no sean estacionarios, por ello usaremos el paquete forecast (Hyndman RJ 2008).

library(forecast)
# prueba para la raíz unitaria y el número de diferencias requeridas, también puede probar la estacionalidad con nsdiffs
ndiffs(chicken, alpha = 0.05, test = c("kpss"))

## [1] 1

ndiffs(egg, alpha = 0.05, test = c("kpss"))

## [1] 1

# series de tiempo diferenciadas
dchick <- diff(chicken)
degg <- diff(egg)
plot.ts(dchick)

plot.ts(degg)

Realizar prueba de causalidad de Granger

Para testear la causalidad en el sentido de granger,la prueba debe correrse varias veces, utilizando un número diferente de retardos en cada ocasión, hasta que aparezca contradicción, es decir que cambie la conclusión de rechazar la hipótesis nula, es decir, puede que con un determinado número de rezagos se rechace la hipótesis nula, pero al aumentar el número de rezagos, no se rechace dicha hipótesis, tambien puede ocurrir que con un número determinado de rezagos la hipotesis nula no se rechace, pero al aumentar el número de rezagos, si se rechace, por eso es importante verificar la prueba en varias ocasiones con un número diferente de rezagos en cada una.

Entonces, en cada prueba al argumento “order” debe especificársele un numero distinto: 1,2,3,4…12, etc.

Para fines didácticos se ha ejecutado la prueba de causalidad en el sentido de Granger con 1,4,6 y 7 rezagos, rechazándose en los primeros tres casos la hipótesis nula y verificándose que desde 7 o más rezagos la hiótesis nula no se rechaza.

library(lmtest)
grangertest(dchick ~ degg, order=1)

## Granger causality test
## 
## Model 1: dchick ~ Lags(dchick, 1:1) + Lags(degg, 1:1)
## Model 2: dchick ~ Lags(dchick, 1:1)
##   Res.Df Df      F   Pr(>F)   
## 1     49                      
## 2     50 -1 10.369 0.002277 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación: Debido a que el p-value es menor que el nivel de significancia, hay evidencia a favor de rechazar la hipótesis nula, por lo tanto podemos concluir que hay evidencia de que la serie de tiempo egg (en el sentido de Granger) causa a la serie de tiempo chick, o en otras palabras quiere decir que conocer los valores de la población de la producción de huevos es valioso para pronosticar los valores futuros de la población de pollos en Estados Unidos. (Zeileis and Hothorn 2002)

library(lmtest)
grangertest(dchick ~ degg, order=4)

## Granger causality test
## 
## Model 1: dchick ~ Lags(dchick, 1:4) + Lags(degg, 1:4)
## Model 2: dchick ~ Lags(dchick, 1:4)
##   Res.Df Df      F   Pr(>F)   
## 1     40                      
## 2     44 -4 4.1762 0.006414 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se verifica que con 4 rezagos, se rechaza la hipótesis nula.

library(lmtest)
grangertest(dchick ~ degg, order=6)

## Granger causality test
## 
## Model 1: dchick ~ Lags(dchick, 1:6) + Lags(degg, 1:6)
## Model 2: dchick ~ Lags(dchick, 1:6)
##   Res.Df Df      F  Pr(>F)  
## 1     34                    
## 2     40 -6 2.6883 0.03034 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Con 6 rezagos se sigue verificando que que se rechaza la hipótesis nula.

library(lmtest)
grangertest(dchick ~ degg, order=7)

## Granger causality test
## 
## Model 1: dchick ~ Lags(dchick, 1:7) + Lags(degg, 1:7)
## Model 2: dchick ~ Lags(dchick, 1:7)
##   Res.Df Df      F  Pr(>F)  
## 1     31                    
## 2     38 -7 2.0445 0.08071 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Con 7 o más rezagos se verifica que No se rechaza la hipótesis nula, debido a que el p-value es mayor al nivel de significancia, es decir que hay evidencia de que la serie de tiempo egg (en el sentido de Granger) NO causa a la serie de tiempo chick.

Es interesante realizar la prueba en sentido inverso para descartar que exista una causalidad bilateral.

Realizar la prueba de causalidad de Granger a la inversa

grangertest(degg ~ dchick, order=1)

## Granger causality test
## 
## Model 1: degg ~ Lags(degg, 1:1) + Lags(dchick, 1:1)
## Model 2: degg ~ Lags(degg, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     49                 
## 2     50 -1 0.5434 0.4646

Interpretación: Debido a que el p-value es mayor que el nivel de significancia, concluimos que existe evidencia a favor de NO rechazar la hipótesis nula, por lo tanto la serie de tiempo de la población de pollos no causa en el sentido de Granger al número de huevos en Estados Unidos.

grangertest(degg ~ dchick, order=4)

## Granger causality test
## 
## Model 1: degg ~ Lags(degg, 1:4) + Lags(dchick, 1:4)
## Model 2: degg ~ Lags(degg, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     40                 
## 2     44 -4 0.2817 0.8881

Se verifica que con 4 rezagos la hipótesis nula no se rechaza.

grangertest(degg ~ dchick, order=6)

## Granger causality test
## 
## Model 1: degg ~ Lags(degg, 1:6) + Lags(dchick, 1:6)
## Model 2: degg ~ Lags(degg, 1:6)
##   Res.Df Df     F Pr(>F)
## 1     34                
## 2     40 -6 0.306 0.9295

Se verifica que con 6 rezagos la hipótesis nula no se rechaza.

grangertest(degg ~ dchick, order=7)

## Granger causality test
## 
## Model 1: degg ~ Lags(degg, 1:7) + Lags(dchick, 1:7)
## Model 2: degg ~ Lags(degg, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     31                 
## 2     38 -7 0.4101 0.8887

Se verifica que con 7 o más rezagos la hipótesis nula NO se sigue rechazando, por lo tanto podemos concluir que en nigún momento hay evidencia de que la serie de tiempo chick (en el sentido de Granger) causa a la serie de tiempo egg.

Por tanto solo existe causalidad unidireccional de la serie de tiempo egg a la serie de tiempo chick, en el sentido de granger.(Esto utilizando como máximo 6 rezagos).

Referencias

Cuellar, Lourdes. 2020. “TEST DE RAIZ UNITARIA DE DICKEY-FULLER EN r.” Youtube. 2020. https://youtu.be/fMqwBJrxJ8s.

Documentation, R. n.d. “Chikens, Eggs and Causality.” http://math.furman.edu/~dcs/courses/math47/R/library/lmtest/html/ChickEgg.html.

Fox, John, and Sanford Weisberg. 2019. An r Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Gujarati, Damodar, and Dawn Porter. 2010. “Econometrı́a (Quinta Edición).” México: Editorial Mc. Graw Hill.

Hyndman RJ, Khandakar. 2008. Automatic Time Series Forecasting: The Forecast Package for r. https://www.jstatsoft.org/article/view/v027i03.

Lesmeister, Cory. 2013. “Chicken or the Egg? Granger-Causality for the Masses.” https://www.r-bloggers.com/2013/06/chicken-or-the-egg-granger-causality-for-the-masses/.

Mendoza González, Miguel Ángel, and Luis Quintana Romero. n.d. “Cointegración y Modelos de Corrección de Error.” https://docs.google.com/viewerng/viewer?url=http%3A%2F%2Fsaree.com.mx%2FeconometriaR%2Fsites%2Fdefault%2Ffiles%2FCap10_MiguelM_LuisQ.pdf&authuser=1.

Pfaff, B. 2008. Analysis of Integrated and Cointegrated Time Series with r. Second. New York: Springer. http://www.pfaffikus.de.

Rico, Victor A. n.d. “Prueba de Raiz Unitaria de Dikey Fuller En r.” https://ricovictor.com/index.php/2021/06/08/prueba-de-raiz-unitaria-de-dickey-fuller-en-r/.

Trapletti, Adrian, and Kurt Hornik. 2021. Tseries: Time Series Analysis and Computational Finance. https://CRAN.R-project.org/package=tseries.

Wickham, Hadley, and Jennifer Bryan. 2019. Readxl: Read Excel Files. https://CRAN.R-project.org/package=readxl.

Wickham, Hadley, Romain Francois, Lionel Henry, and Kirill Muller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Zeileis, Achim, and Torsten Hothorn. 2002. “Diagnostic Checking in Regression Relationships.” R News 2 (3): 7–10. https://CRAN.R-project.org/doc/Rnews/.

Zhu, Hao. 2021. kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.