AR(1): Autoregressive porcess of order 1

Un prioceso es AR(1), si se puede escribir como:

\[\begin{align*} y_{t}&=\phi y_{t-1} + \epsilon_{t} \quad \epsilon_{t}\overset{iid}{\sim} N(0,v)\\ y_{t-1}&=\phi y_{t-2} + \epsilon_{t-1}, \quad \text{de manera recursiva tenemos que}\\ &=\phi^ky_{t-k}+\sum_{j=0}^{k-1}\phi^j\epsilon_{t-j} \end{align*}\] Con \(\phi\) el coeficiente de autoregresión.

Sí \(\phi \in (-1,1)\) entonces \[ y_t= \sum_{j=0}^{\infty}\phi^j\hspace{0.1cm}\epsilon_{t-j} \]

Propiedades del AR(1)

  1. \[E(y_t)=0 \quad\text{y} \quad Var(y_t)= \sum_{j=0}^{\infty} \phi^{2j}v=\frac{v}{1-\phi^2}\]

  2. Si \(\phi \in (-1,1)\) el proceso va a ser estacionario.

Función de autocovarianza

  1. \(\gamma(h)= \mathbf{E}[y_t \hspace{0.1cm} y_{t-h}]=\mathbf{E}[\sum_{j=0}^{\infty}\phi^j\hspace{0.1cm}\epsilon_{t-j} \cdot \sum_{k=0}^{\infty}\phi^k\hspace{0.1cm}\epsilon_{t-h-k} ] =v\phi^h\sum_{j=0}^{\infty}\phi^{2j}\)

De nuevo veamos que sí \(\phi \in (-1,1)\) tenemos qeu \[\gamma(h)= \frac{v\phi^h}{1-\phi^2}\]

Función de autocorrelación: PACF Partial Autocorrelation Function

  1. \[ \rho(h)=\frac{\gamma(h)}{\gamma(0)}\] sí \(\phi \in (-1,1)\) se tiene que \(\rho(h)=\phi^h\) De manera general \(\forall h \in \mathbb{Z}\)

\[ \rho(h)=\phi^{\lvert{h}\rvert} \qquad \gamma(h)=\frac{v\phi^{\lvert {h}\rvert}}{1-\phi^2} \]

Simulando un proceso AR(1)

###############################
#####    sample ar(1)     #####
###############################
#
# sample data from 2 ar(1) processes and plot their ACF and PACF functions
#
set.seed(2021)
T=500 # number of time points
#
# sample data from an ar(1) with ar coefficient phi = 0.9 and variance 1
#
v=1.0 # innovation variance
sd=sqrt(v) #innovation stantard deviation
phi1=0.9 # ar coefficient

#ARIMA: Autoregressive integrated moving average process:
yt1=arima.sim(n = T, model = list(ar = phi1), sd = sd)
#como no tenemos parte de moovign average ni un proceso       integradomentonces solo ponemos en formato lista.

Ahora vamos a simular un AR(1), cuyo \(\phi\) es el inverso aditivo del otro es decir \(\phi=-0.9\) y graficarlos

# sample data from an ar(1) with ar coefficient phi = -0.9 and variance 1
#
phi2=-0.9 # ar coefficient
yt2=arima.sim(n = T, model = list(ar = phi2), sd = sd)

par(mfrow = c(2, 1), cex.lab = 1.3)
plot(yt1,main=expression(phi==0.9))
plot(yt2,main=expression(phi==-0.9))

Ahora vamos a graficar la función de autocovariabza ACF real y la muestral:

ACF real:

par(mfrow = c(1, 2), cex.lab = 1.3)
lag.max=50 # max lag
#
## plot true ACFs for both processes
#
cov_0=sd^2/(1-phi1^2) # compute auto-covariance at h=0
cov_h=phi1^(0:lag.max)*cov_0 # compute auto-covariance at h
plot(0:lag.max, cov_h/cov_0, pch = 1, type = 'h', col = 'red',
     ylab = "true ACF", xlab = "Lag",ylim=c(-1,1), main=expression(phi==0.9))

cov_0=sd^2/(1-phi2^2) # compute auto-covariance at h=0
cov_h=phi2^(0:lag.max)*cov_0 # compute auto-covariance at h
# Plot autocorrelation function (ACF)
plot(0:lag.max, cov_h/cov_0, pch = 1, type = 'h', col = 'red',
     ylab = "true ACF", xlab = "Lag",ylim=c(-1,1),main=expression(phi==-0.9))

Función de autocovarianza muestral:

par(mfrow = c(1, 2), cex.lab = 1.3)
## plot sample ACFs for both processes
#
acf(yt1, lag.max = lag.max, type = "correlation", ylab = "sample ACF",
    lty = 1, ylim = c(-1, 1), main = " ")
acf(yt2, lag.max = lag.max, type = "correlation", ylab = "sample ACF",
    lty = 1, ylim = c(-1, 1), main = " ")

Función de autocorrelación parcial (PACF) muestral:

par(mfrow = c(1, 2), cex.lab = 1.3)
## plot sample PACFs for both processes
#
pacf(yt1, lag.ma = lag.max, ylab = "sample PACF", ylim=c(-1,1),main="")
pacf(yt2, lag.ma = lag.max, ylab = "sample PACF", ylim=c(-1,1),main="")

Sabemos que \(\texttt{PACF}=\phi\) para \(n=1\) y que para \(n>=2\) la\(\texttt{PACF}=0\) , entonces deberiamos ver que todas las barras se acerquen a 0, a medida que \(n\) aumente.

NOTAS:

  1. LA ACF de un AR(1) con coeficiente \(\phi>0\) decae exponencialmente.
  2. LA ACF de un AR(1) con coeficiente \(\phi<0\) decae exponencialmente de forma oscilatoria
  3. Los coeficientes PACF para rezagos superiores a 1 (n>1) son cero
  4. El coeficiente PACF en el retardo 1 (n=1), \(\phi(1,1)\) es igual a \(\phi\)

Ejercicios

  1. ¿Cuál de las siguientes corresponde a la función de autocovarianza en el retardo \(h=2, \gamma(2)\) del proceso autorregresivo: \(Y_t=0.7Y_{t-1}+\epsilon_t, \quad \epsilon_t \overset{iid}{\sim} N(0,v)\), con \(v=2\)

Como \(\phi=0.7>0\) tenemos que

\[ \gamma(h=2)=\frac{v\phi^{|h|}}{1-\phi^2}=\frac{2(0.7)^{|2|}}{1-0.7^2}=1.9216 \]

  1. ¿Cuál es el coeficiente PACF en el retardo 1 para el proceso AR(1) \(Y_t=-0.7Y_{t-1}+\epsilon_t, \quad \epsilon_t \overset{iid}{\sim} N(0,1)\)?

Como el retardo(lag) es 1, es decir \(n=1\), tenemos que \(PACF=\phi=-0.7\)

  1. ¿Cuál es la función de autovarianza en el retardo 1, $ (1)$) del proceso AR(1) \(y_t=0.6Y_{t-1}+\epsilon_t, \quad \epsilon_t \overset{iid}{\sim} N(0,v)\) con varianza \(v=2\)

\[ \gamma(h=1)=\frac{v\phi^{|h|}}{1-\phi^2}=\frac{2(0.6)^{|1|}}{1-0.6^2}=1.875 \]

MODELOS DE REGRESIÓN

Estimación de máxima verosimilitud

sea \(\mathbf{y}=\mathbf{X}\mathbf{\beta}+\mathbf{\epsilon}\), $N(,v) $ \[ \mathbf{\beta}_{MLE}=(\mathbf{X^TX})^{-1}\mathbf{X^Ty } \] con un estimador insesgado para la varianza: \[ \hat{v}=s^2=\frac{(\mathbf{y}-\mathbf{X\hat{\beta}})^T(\mathbf{y}-\mathbf{X\hat{\beta}})}{n-p} \]

INFERENCIA BAYESIANA

La función de verosimilitud es:

\[ p(y | \beta, v) = \frac{1}{(2\pi v)^{n/2}} \exp\left\{-\frac{1}{2} (y - X\beta)^T (y - X\beta)\right\} \]

Si se utiliza una prior de la forma \(p(\beta, v) \propto \frac{1}{v}\), la distribución posterior se expresa como: \[\begin{equation} p(\beta, v | y) \propto \frac{1}{v^{n/2 + 1}} \exp\left\{-\frac{1}{2v} (y - X\beta)^T (y - X\beta)\right\} \end{equation}\]

Además, se puede demostrar que: \[\begin{align} p(\beta | v, y) &\sim \mathcal{N}(\tilde{\beta}_{\text{MLE}}, v (X^TX)^{-1}) \\ p(v | y) &\sim \text{IG}\left(\frac{n - p}{2}, \frac{d}{2}\right) \end{align}\] donde \(d = (y - X\hat{\beta}_{\text{MLE}})^T (y - X\hat{\beta}_{\text{MLE}})\).

con \(p=dim(\beta)\)

Dado que \(p(\beta, v | y)=p(\beta | v, y)p(v | y)\) podemos muestrear a partir de la distribución posterior de \(\beta\) y \(v\)

Estimación de máxima verosimilitud en el AR(1)

Sea un proceso AR(1): \[ y_{t}=\phi y_{t-1} + \epsilon_{t} \quad \epsilon_{t}\overset{iid}{\sim} N(0,v) \quad \phi \in(-1,1) \] Entonces el proceso es estacionario y todas las \(y_t\) son estacionarias: tenemos que parala primera es: \(y_1 \sim N(0,\frac{v}{1-\phi^2})\) y para el resto será: \[y_t|y_{t-1}\sim N(\phi y_{t-1},v)\] Donde tenemos que su distribución se puede escribir para 1:T, como sigue \[ \begin{align} p(y_{1:T}|\phi,v) &= p(y_{1}|\phi,v) \prod_{t=2}^{T} p(y_t|y_{t-1},\phi,v) \\ &= \frac{(1-\phi^2)^{1/2}}{(2\pi v)^{T/2}} \exp\left\{-\frac{Q^{*}(\phi)}{2v} \right\} \end{align} \] Donde \(Q^{*}(\phi)=y^2_1(1-\phi^2) + \sum_{t=2}^T(y_t-\phi y_{t-1})^2\) Y vamos a denotar \(Q(\phi)=\sum_{t=2}^T(y_t-\phi y_{t-1})^2\)

Aunque conocer la distribución de 1:T , no nos es tan útil dado que al conocer la distribución de \(y_1\), es mejor obtener la distribución de 2:T , usando la verosimilitud condicional

De manera gráfica tenemos que, la regresión para un proceso AR(1) es: \[ \begin{align} p(y_{1:T}|\phi,v) &= \frac{1}{(2\pi v)^{T- 1/2}} \exp\left\{-\frac{Q(\phi)}{2v} \right\} \end{align} \] \[ \begin{bmatrix} y_{2} \\ \vdots \\ y_{T} \end{bmatrix} = \begin{bmatrix} y_{1} \\ \vdots \\ y_{T-1} \\ \end{bmatrix} \phi + \begin{bmatrix} \epsilon_{1} \\ \vdots \\ \epsilon_{T} \\ \end{bmatrix} \\ y=\mathbf{X}\mathbf{\beta}+\mathbf{\epsilon} \quad N(0,v\mathbf{I}) \]

Que se asemeja mucho a un modelo de regresión, donde \(\phi\) es \(\beta\) entonces para obtener un estimador \(\texttt{MLE}\) de \(\phi\) tenemos que: \[ \hat{\phi}_{MLE}= \frac{\sum_{t=2}^T(y_t\cdot y_{t-1})}{\sum_{t=2}^T(y^2_{t-1})} \\ \hat{v}=s^2=\frac{\sum_{t=2}^T(y_t-\hat{\phi}_{MLE} \cdot y_{t-1})^2}{T-2} \]

En el caso de la verosimilitud completa no es posible establecer una correspondencia entre las ecuaciones que tenemos para la verosimilitud del proceso autorregresivo y el modelo normal lineal, entonces hay que hacer optimización numerica para obtener \(\phi\) y \(v\)

Entonces, como no podemos usar los mismos pasos de la regesión normal lineal para econtrar el \(\phi\) que maximize la función de verosimilitud: \(p(y_{1:T}|\phi,v)\), entonces voy a usar el algoritmo de Newton Rhapson para maximizar $()= $ y obtener el estimador maximo verosimil de \(\phi\)

El siguiente código le permite calcular el MLE del coeficiente AR \(\phi\), el estimador insesgado de \(\sigma^2\), \(\sigma^2_s\), y el MLE de \(\sigma^2\) basándose en un conjunto de datos simulado a partir de un proceso AR(1) y utilizando la verosimilitud condicional.

####################################################
#####             MLE for AR(1)               ######
####################################################
set.seed(2021)
phi=0.9 # ar coefficient
v=1
sd=sqrt(v) # innovation standard deviation
T=500 # number of time points
yt=arima.sim(n = T, model = list(ar = phi), sd = sd) 

## Case 1: Conditional likelihood
y=as.matrix(yt[2:T]) # response
X=as.matrix(yt[1:(T-1)]) # design matrix
phi_MLE=as.numeric((t(X)%*%y)/sum(X^2)) # MLE for phi
s2=sum((y - phi_MLE*X)^2)/(length(y) - 1) # Unbiased estimate for v 
v_MLE=s2*(length(y)-1)/(length(y)) # MLE for v

cat("\n MLE of conditional likelihood for phi: ", phi_MLE, "\n",
    "MLE for the variance v: ", v_MLE, "\n", 
    "Estimate s2 for the variance v: ", s2, "\n")
## 
##  MLE of conditional likelihood for phi:  0.9261423 
##  MLE for the variance v:  1.048 
##  Estimate s2 for the variance v:  1.050104

Este código le permite calcular estimaciones del coeficiente AR(1)\(\phi\) y la varianza utilizando la función arima en R. El primer caso utiliza la suma condicional de cuadrados; el segundo y el tercero utilizan la verosimilitud total con diferentes puntos de partida para la optimización numérica necesaria para calcular el MLE con la verosimilitud t

# Obtaining parameter estimates using the arima function in R

set.seed(2021)
phi=0.9 # ar coefficient
v=1
sd=sqrt(v) # innovation standard deviation
T=500 # number of time points
yt=arima.sim(n = T, model = list(ar = phi), sd = sd) 

#Using conditional sum of squares, equivalent to conditional likelihood 
arima_CSS=arima(yt,order=c(1,0,0),method="CSS",n.cond=1,include.mean=FALSE)
cat("AR estimates with conditional sum of squares (CSS) for phi and v:", arima_CSS$coef,arima_CSS$sigma2,
"\n")
## AR estimates with conditional sum of squares (CSS) for phi and v: 0.9261423 1.048
#Uses ML with full likelihood 
arima_ML=arima(yt,order=c(1,0,0),method="ML",include.mean=FALSE)
cat("AR estimates with full likelihood for phi and v:", arima_ML$coef,arima_ML$sigma2,
"\n")
## AR estimates with full likelihood for phi and v: 0.9265251 1.048434
#Default: uses conditional sum of squares to find the starting point for ML and 
#         then uses ML 
arima_CSS_ML=arima(yt,order=c(1,0,0),method="CSS-ML",n.cond=1,include.mean=FALSE)
cat("AR estimates with CSS to find starting point for ML for phi and v:", 
arima_CSS_ML$coef,arima_CSS_ML$sigma2,"\n")
## AR estimates with CSS to find starting point for ML for phi and v: 0.9265252 1.048434

Este código le muestra cómo calcular la MLE para \(\phi\) utilizando la verosimilitud total y la función optimize en R

set.seed(2021)
phi=0.9 # ar coefficient
v=1
sd=sqrt(v) # innovation standard deviation
T=500 # number of time points
yt=arima.sim(n = T, model = list(ar = phi), sd = sd) 

## MLE, full likelihood AR(1) with v=1 assumed known 
# log likelihood function
log_p <- function(phi, yt){
  0.5*(log(1-phi^2) - sum((yt[2:T] - phi*yt[1:(T-1)])^2) - yt[1]^2*(1-phi^2))
}

# Use a built-in optimization method to obtain maximum likelihood estimates
result =optimize(log_p, c(-1, 1), tol = 0.0001, maximum = TRUE, yt = yt)
cat("\n MLE of full likelihood for phi: ", result$maximum)
## 
##  MLE of full likelihood for phi:  0.9265928

Estimación de bayesiana en el AR(1)