Autoregresivo

2023

Ejemplo 1: Longley

Para ilustrar la metodología usaremos una base de datos incorporada en R llamada datos de la regresión de Longley donde la variable respuesta es el número de personas empleadas anualmente (Employed) de 1947 hasta 1962 y utilizaremos como variables predictoras el deflactor del precio implícito (GNP y en español es el PIB) con índice 1954=100 y la población (Population) no institucionalizada de 14 o más años de edad. Los datos aparecieron originalmente en Longley (1967)[https://www.jstor.org/stable/2283673].

Ejemplo 1: Longley

Uno de los principales problemas en la economía es evaluar el crecimiento económico, es decir, el valor de los bienes y servicios producidos por una economía a lo largo del tiempo. El principal escollo para medir esto es la distorsión que genera el incremento de los precios y su repercusión sobre el valor de lo producido, conocido como inflación.

Ejemplo 1: Longley

Para evaluar el crecimiento real, y no sólo de su valor, se utiliza el producto interior bruto real, en el que se consideran sólo las variaciones de las cantidades producidas en términos reales. Para ello se hace necesario eliminar el efecto de los precios sobre la economía. Por tanto es necesario el ajuste a través de un deflactor.
Un deflactor es un índice de precios, simple o compuesto, que permite desagregar las series en sus dos componentes de precios y cantidades.

Ejemplo 1: Longley

El ejemplo más útil es el deflactor del PIB que resulta ser el cociente entre el PIB nominal y el PIB real, expresado en forma de Índice, o la relación entre el valor del PIB del año en curso y el valor del PIB del año base. El deflactor del PIB entonces sirve para medir el crecimiento económico real -no meramenente nominal- de una economía.

Ejemplo 1: Longley

longley
##      GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
## 1947         83.0 234.289      235.6        159.0    107.608 1947   60.323
## 1948         88.5 259.426      232.5        145.6    108.632 1948   61.122
## 1949         88.2 258.054      368.2        161.6    109.773 1949   60.171
## 1950         89.5 284.599      335.1        165.0    110.929 1950   61.187
## 1951         96.2 328.975      209.9        309.9    112.075 1951   63.221
## 1952         98.1 346.999      193.2        359.4    113.270 1952   63.639
## 1953         99.0 365.385      187.0        354.7    115.094 1953   64.989
## 1954        100.0 363.112      357.8        335.0    116.219 1954   63.761
## 1955        101.2 397.469      290.4        304.8    117.388 1955   66.019
## 1956        104.6 419.180      282.2        285.7    118.734 1956   67.857
## 1957        108.4 442.769      293.6        279.8    120.445 1957   68.169
## 1958        110.8 444.546      468.1        263.7    121.950 1958   66.513
## 1959        112.6 482.704      381.3        255.2    123.366 1959   68.655
## 1960        114.2 502.601      393.1        251.4    125.368 1960   69.564
## 1961        115.7 518.173      480.6        257.2    127.852 1961   69.331
## 1962        116.9 554.894      400.7        282.7    130.081 1962   70.551

Ejemplo 1: Longley

El modelo lineal ajustado fue el siguiente:

data(longley) 
g <- lm(Employed ~ GNP + Population, data=longley)
summary(g,cor=T)
## 
## Call:
## lm(formula = Employed ~ GNP + Population, data = longley)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.80899 -0.33282 -0.02329  0.25895  1.08800 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 88.93880   13.78503   6.452 2.16e-05 ***
## GNP          0.06317    0.01065   5.933 4.96e-05 ***
## Population  -0.40974    0.15214  -2.693   0.0184 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5459 on 13 degrees of freedom
## Multiple R-squared:  0.9791, Adjusted R-squared:  0.9758 
## F-statistic: 303.9 on 2 and 13 DF,  p-value: 1.221e-11
## 
## Correlation of Coefficients:
##            (Intercept) GNP  
## GNP         0.98            
## Population -1.00       -0.99
##             (Intercept)           GNP Population
## (Intercept) 190.0269691  0.1445617813 -2.0954381
## GNP           0.1445618  0.0001133631 -0.0016054
## Population   -2.0954381 -0.0016053999  0.0231456
##             (Intercept)        GNP Population
## (Intercept)   1.0000000  0.9849405 -0.9991552
## GNP           0.9849405  1.0000000 -0.9910901
## Population   -0.9991552 -0.9910901  1.0000000

Ejemplo 1: Longley

Comparando la correlación entre las variables GNP y Population y sus correpondientes coeficientes. Podemos apreciar que están altamente correlacionados negativamente.
En datos colectados a través del tiempo tales como estos, los errores sucesivos podrían estar correlacionados. Asumiendo que los errores toman una forma autoregresiva simple: \[\varepsilon_{i+1}=\rho \varepsilon_{i}+\delta_i \quad \delta_i \sim N(0,\tau^2)\]

Ejemplo 1: Longley

Ejemplo 1: Longley

Podemos estimar esta correlación \(\rho\) por

cor(g$res[-1],g$res[-16])
## [1] 0.3104092

El grafico de los residuales a través del tiempo es el siguiente:

res<-g$res
plot(seq(1:16),res, pch=16, xlab="t", ylab=expression(widehat(paste(epsilon))))
abline(h=0,lty=2)

Ejemplo 1: Longley

Ejemplo 1: Longley

Por simplicidad, se puede construir la matriz de covarianza del error \(\textbf{V}_{ij}=\rho^{|i-j|}\). Asumiendo que se conoce \(\rho=0.310409\). Luego \(\textbf{V}\) es calculada de la siguiente forma: 

V <- diag(16)
V <- 0.310409^abs(row(V)-col(V))
V

Ejemplo 1: Longley

Ejemplo 1: Longley

y el estimador de mínimos cuadrados generalizados para \(\beta\) es dado por \(\widehat{\boldsymbol{\beta}}=\left(\boldsymbol{X}^{\top} \textbf{V}^{-1} \boldsymbol{X}\right)^{-1}\boldsymbol{X}^{\top} \textbf{V}^{-1} \boldsymbol{y}\) es 

g     <- lm(Employed ~ GNP + Population, data=longley)
X     <- model.matrix(g)
V.inv <- solve(V)
beta  <- solve(t(X)%*%V.inv%*%X)%*%t(X)%*%V.inv%*%longley$Empl
beta

Ejemplo 1: Longley

El error estandar de \(\widehat{\boldsymbol{\beta}}\), \(\sqrt{Var(\widehat{\boldsymbol{\beta}})}=\sqrt{\widehat{\sigma}^2 \left( \boldsymbol{X}^{\top} \textbf{V}^{-1} \boldsymbol{X}\right)^{-1}}\) es 

sigma_est <- (t(longley$Empl)%*%V.inv%*%longley$Empl-t(beta)%*%t(X)%*%V.inv%*%longley$Empl)/g$df
sqrt(diag(solve(t(X)%*%V.inv%*%X))*sigma_est)
## (Intercept)         GNP  Population 
## 19.07759167  0.01464312  0.20984390

Ejemplo 1: Longley

Otra opción es usar el paquete nlme(), que contiene la función gls(). Podemos usar dicha función de la siguiente forma: 

library(nlme)
g1<-gls(Employed~GNP+Population,correlation=corAR1(form=~Year),data=longley)
summary(g1)

Ejemplo 1: Longley

Ejemplo 1: Longley

Se observa que el valor estimado de \(\rho\) es 0.64. Sin embargo, al encontrar los intervalos de confianza para \(\rho\): 

intervals(g1)

Ejemplo 1: Longley

Prueba de hipotesis de Durbin-Watson

Para realizar la prueba de hipótesis de Durbin Watson necesitamos extraer los residuales del modelo calculado usando mínimos cuadrados ordinarios.

\[H_0: \rho=0\] \[H_1: \rho>0\] o

\[H_0: DW=2\] \[H_0: DW<2\]

Prueba de hipotesis de Durbin-Watson

La estadística de prueba es dada por:

\[DW=\frac{\sum_{i=2}^{n} (\widehat{e}_t-\widehat{e}_{t-1})^2}{\sum_{i=1}^{n} \widehat{e}^2_t}\]

g <- lm(Employed ~ GNP + Population, data=longley)
residuales <- g$res
residuales
##        1947        1948        1949        1950        1951        1952 
##  0.67521130  0.30582253 -0.09098821 -0.27823770 -0.57801231 -0.80898950 
##        1953        1954        1955        1956        1957        1958 
##  0.12689319 -0.49655508  0.07001902  1.08799625  0.61089180 -0.54070252 
##        1959        1960        1961        1962 
## -0.22904034  0.24332304  0.04440329 -0.14203475

Prueba de hipotesis de Durbin-Watson

Prueba de hipotesis de Durbin-Watson

Los cálculos de la tabla puden ser obtenidos utilizando el siguiente código de R: 

g <- lm(Employed ~ GNP + Population, data=longley)
dif <- c()
dif[1] <- 0
for(i in 2:length(longley$Empl))
{
  dif[i] <- g$res[i]-g$res[i-1]
}

DWtest <- sum(dif^2)/sum(g$res^2)
DWtest
## [1] 1.301484

Prueba de hipotesis de Durbin-Watson

Luego usando las tablas estadísticas, donde \(k\) es el nùmero de regresores exluyendo el intercepto:

Prueba de hipotesis de Durbin-Watson

Analizando los resultado tenemos que:

Prueba de hipotesis de Durbin-Watson

Luego usando el paquete lmtest 

library(lmtest)
g <- lm(Employed ~ GNP + Population, data=longley)
dwtest(g)
## 
##  Durbin-Watson test
## 
## data:  g
## DW = 1.3015, p-value = 0.02245
## alternative hypothesis: true autocorrelation is greater than 0

Como el valor de DW obtenido a partir de los residuos (1.3015) está en la región de indecisión, la prueba sería inconclusa, es decir, no hay evidencias, para afirmar que los errores sean o no autocorrelacionados.

Prueba de hipotesis de Q de Ljung-Box

\(H_{0}\) :las autocorrelaciones hasta el retraso \(k\) son todas 0
\(H_{1}\) :las autocorrelaciones de uno o más retrasos difieren de 0

Box.test(residuals(g), lag = 1, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(g)
## X-squared = 1.5905, df = 1, p-value = 0.2073

Ejemplo 2: Hartnagel

  • Examinemos los datos de series temporales sobre las tasas de delincuencia femenina en Canadá analizados por Fox y Hartnagel (1979). El conjunto de datos se llama Hartnagel y se encuentra en el paquete carData pero que al cargar el paquete car se carga también:
library("car")
## Loading required package: carData
data(Hartnagel)

Ejemplo 2: Hartnagel

  • Las varibles del conjunto de datos son:

    • year: 1931 - 1968
    • tfr: la tasa de fertilidad total, nacimientos por 1000 mujeres
    • partic: tasa de participación de las mujeres en la fuerza laboral, por 1000
    • degrees: tasa de estudios postsecundarios de mujeres, por 100000
    • fconvict: tasa de condena por delito procesable de las mujeres, por 100000

Ejemplo 2: Hartnagel

  • ftheft: tasa de condena por robo de mujeres, por 100000
  • mconvict: tasa de condena por delito procesable de los hombres, por 100000.
  • mtheft: tasa de condena por robo de hombres, por 100000

Ejemplo 2: Hartnagel

  • Ajustaremos la regresión considerando a fconvict como la variable dependiente y tfr, partic, degrees y mconvictcomo los regresones.

  • La razón para incluir el último predictor es controlar las variables omitidas que afectan la tasa de criminalidad en general.

Ejemplo 2: Gráfico

plot(fconvict ~ year, type="n",data=Hartnagel,
ylab="Convictions per 100,000 Women")
grid(lty=1)
with(Hartnagel, points(year, fconvict, type="o", pch=16))

Ejemplo 2: Gráfico

Ejemplo 2: Ajuste

  • Podemos ver que la tasa de convicción de las mujeres fluctuó sustancialmente pero gradualmente durante este período histórico, sin una tendencia general aparente.
  • Ajustamos una regresión por medio del métodod de MCO que produce el siguiente resultado:
mod.ols <- lm(fconvict ~ tfr + partic + degrees + mconvict, data=Hartnagel)
summary(mod.ols)

Ejemplo 2: Ajuste

## 
## Call:
## lm(formula = fconvict ~ tfr + partic + degrees + mconvict, data = Hartnagel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.964  -9.204  -3.566   6.149  48.385 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 127.639997  59.957044   2.129   0.0408 *  
## tfr          -0.046567   0.008033  -5.797 1.75e-06 ***
## partic        0.253416   0.115132   2.201   0.0348 *  
## degrees      -0.212049   0.211454  -1.003   0.3232    
## mconvict      0.059105   0.045145   1.309   0.1995    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.19 on 33 degrees of freedom
## Multiple R-squared:  0.6948, Adjusted R-squared:  0.6578 
## F-statistic: 18.78 on 4 and 33 DF,  p-value: 3.905e-08

Ejemplo 2: Análisis de la salida

  • La tasa de criminalidad de las mujeres, por lo tanto, parece disminuir con la fertilidad y aumentar con la participación de la fuerza laboral.

  • Los coeficientes para los otros dos predictores tienen p-values grandes.

Ejemplo 2: Gráfico

  • Sin embargo, un gráfico de los residuos de la regresión de MCO sugiere que pueden estar sustancialmente autocorrelacionados:
plot(Hartnagel$year, residuals(mod.ols), type="o", pch=16,
xlab="Year", ylab="OLS Residuals")
abline(h=0, lty=2)

Ejemplo 2: Gráfico

Ejemplo 2: Hartnagel

  • La función acf() en el paquete de R stats calcula y grafoca las funciones de autocorrelación y autocorrelación parcial de una serie de tiempo, aquí para los residuos obtenidos por MCO, porque los residuos varían sistemáticamente con el tiempo:
acf(residuals(mod.ols))
acf(residuals(mod.ols), type="partial")

Ejemplo 2: Hartnagel

Ejemplo 2: Hartnagel

Ejemplo 2: Hartnagel

  • Las líneas horizontales discontinuas (de color azul) en los plots corresponden a límites de confianza aproximados del 95%.

  • El patrón general de las funciones de autocorrelación y autocorrelación parcial: descomposición sinusoidal en la primera y se observan dos picos, uno positivo (\(\phi_1 > 0\)), el otro negativo (\(\phi_2 < 0\)), en el último, sugiere un proceso AR(2).

\[\varepsilon_t = \phi_1 \varepsilon_{t-1} + \phi_2 \varepsilon_{t-2} + v_t\]

Ejemplo 2: Hartnagel

  • Calculando las estadísticas de Durbin-Watson, usando la función durbinWatsonTest() en el paquete car.

  • Por defecto, esta función calcula los p-values por boostrap para las estadísticas de Durbin-Watson.

durbinWatsonTest(mod.ols, max.lag=5)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6883450     0.6168636   0.000
##    2       0.1922665     1.5993563   0.146
##    3      -0.1685699     2.3187448   0.312
##    4      -0.3652775     2.6990538   0.008
##    5      -0.3673240     2.6521103   0.016
##  Alternative hypothesis: rho[lag] != 0

Ejemplo 2: Hartnagel

  • Tres de las cinco estadísticas de Durbin-Watson tienen p-values pequeñas, incluyendo el primero.

  • Otra alternativa, la función dwtest() en el paquete lmtest calcula el p-valor para el estadístico de Durbin-Watson de primer orden y como observamos la prueba de DW es significativo.

library("lmtest")
dwtest(mod.ols, alternative="two.sided")

Ejemplo 2: Hartnagel

## 
##  Durbin-Watson test
## 
## data:  mod.ols
## DW = 0.61686, p-value = 1.392e-08
## alternative hypothesis: true autocorrelation is not 0

Ejemplo 2: Hartnagel

  • Varios de los argumentos para la función gls() son los mismos que los de lm(), en particular, gls() toma argumentos como model, data, subset.

    • El argumento weights en gls() puede ser usado para especificar un modelo para el error de la varianza.
    • El argumento correlation puede ser usado para especificar un modelo para la autocorrelación del error.
    • El argumento method selecciona el método de la estimación method="ML" para la estimación de máxima verosimilitud.

Ejemplo 2: Hartnagel

library("nlme")
mod.gls <- gls(fconvict ~ tfr + partic + degrees + mconvict, data=Hartnagel, correlation=corARMA(p=2), method="ML")
summary(mod.gls)

Ejemplo 2: Hartnagel

## Generalized least squares fit by maximum likelihood
##   Model: fconvict ~ tfr + partic + degrees + mconvict 
##   Data: Hartnagel 
##        AIC      BIC    logLik
##   305.4145 318.5152 -144.7073
## 
## Correlation Structure: ARMA(2,0)
##  Formula: ~1 
##  Parameter estimate(s):
##       Phi1       Phi2 
##  1.0683473 -0.5507269 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 83.34028  59.47084  1.401364  0.1704
## tfr         -0.03999   0.00928 -4.308632  0.0001
## partic       0.28761   0.11201  2.567653  0.0150
## degrees     -0.20984   0.20658 -1.015757  0.3171
## mconvict     0.07569   0.03501  2.161899  0.0380
## 
##  Correlation: 
##          (Intr) tfr    partic degres
## tfr      -0.773                     
## partic   -0.570  0.176              
## degrees   0.093  0.033 -0.476       
## mconvict -0.689  0.365  0.047  0.082
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4991516 -0.3716988 -0.1494540  0.3372409  2.9094711 
## 
## Residual standard error: 17.70228 
## Degrees of freedom: 38 total; 33 residual

Ejemplo 2: Hartnagel

  • Especificando la estructura de correlación como correlation=corARMA(p=2) ajusta un proceso AR(s) para los errores. Noten que aquí el componente de media movil no es considerado pues su orden \(q=0\) y por tanto es ausente.

Ejemplo 2: Hartnagel

  • En este caso, las estimaciones de ML de los parámetros de regresión bajo el modelo de correlación de errores AR(2) no son terriblemente diferentes de las estimaciones de MCO, aunque el coeficiente para mconvict ahora tiene un p-valor más pequeño.

Ejemplo 2: Hartnagel

  • El estimado de ML de los parámetros de los errores autoregresivos son considerables: \(\hat{\phi}_1=1.068\) y \(\hat{\phi}_2=-0.551\).

  • Podemos emplear pruebas de razón de verosimilitud para verificar si los parámetros del proceso AR(2) para los errores son necesarios, y si un modelo autorregresivo de segundo orden es suficiente.

Ejemplo 2: Hartnagel

  • Procedemos actualizando el modelo original ajustado con el gls, reespecificando el proceso de series de tiempo en los errores, luego comparamos los modelos anidados usando la función anova().
mod.gls.3 <- update(mod.gls, correlation=corARMA(p=3))
mod.gls.1 <- update(mod.gls, correlation=corARMA(p=1))
mod.gls.0 <- update(mod.gls, correlation=NULL)

Ejemplo 2: Hartnagel

anova(mod.gls, mod.gls.1) # AR(2) vs AR(1)
##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## mod.gls       1  8 305.4145 318.5152 -144.7073                        
## mod.gls.1     2  7 312.4234 323.8865 -149.2117 1 vs 2 9.008881  0.0027

Ejemplo 2: Hartnagel

anova(mod.gls, mod.gls.0) # AR(2) vs errores no correlacionados
##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## mod.gls       1  8 305.4145 318.5152 -144.7073                        
## mod.gls.0     2  6 339.0011 348.8266 -163.5006 1 vs 2 37.58658  <.0001

Ejemplo 2: Hartnagel

anova(mod.gls.3, mod.gls) # AR(3) vs AR(2)
##           Model df      AIC      BIC    logLik   Test    L.Ratio p-value
## mod.gls.3     1  9 307.3961 322.1343 -144.6980                          
## mod.gls       2  8 305.4145 318.5152 -144.7073 1 vs 2 0.01846713  0.8919
  • En conclusión nos inclinamos con el AR(2) para los errores.