AR(2) model

\(h\)-steps ahead forecasting

1 Data and model

Let us assume that our data set is \(\{y_1,\ldots,y_n\}\) and is observed chronologically. We want to entertain an AR(2) model for the data, ie. \[ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \epsilon_t \qquad \epsilon_t \ \ iid \ \ N(0,\sigma^2), \] where \(\theta=(\phi_1,\phi_2,\sigma)\).

2 Simulating some data

set.seed(1233)
n    = 300
sig  = 1
phi1 = 1.75;phi2 = -0.80
#phi1 = 0.90;phi2 = 0.05
y    = rep(0,2*n) 
for (t in 3:(2*n)) 
  y[t] = phi1*y[t-1]+phi2*y[t-2]+rnorm(1,0,sig) 
y = y[(n+1):(2*n)]

ts.plot(y)

3 h-steps ahead forecasting distributions

Here we are interested in obtaining the forecasting distribution \(h\) steps ahead,i.e. \[ p(y_{n+h}|y_{1:n},\theta) \equiv p_N(y_{n+h}|\mu_{n+h},\sigma_{n+h}^2), \qquad h=1,\ldots,H. \] For \(h=1\), it is easy to see that \[\begin{eqnarray*} y_{n+1} &=& \phi_1 y_n + \phi_2 y_{n-1} + \epsilon_{n+1}\\ E(y_{n+1}|y_{1:n},\theta) &=& \phi_1 y_n + \phi_2 y_{n-1}\\ V(y_{n+1}|y_{1:n},\theta) &=& V(\epsilon_{n+1})=\sigma^2 \end{eqnarray*}\]

For \(h=2\), we first write \(y_{n+2}\) as a function of \((y_n,y_{n-1})\): \[\begin{eqnarray*} y_{n+2} &=& \phi_1 y_{n+1} + \phi_2 y_n + \epsilon_{n+2}\\ &=& \phi_1\{\phi_1 y_n + \phi_2 y_{n-1} + \epsilon_{n+1}\} + \phi_2 y_n + \epsilon_{n+2}\\ &=& (\phi_1^2+\phi_2)y_n + \phi_1\phi_2 y_{n-1} + \phi_1 \epsilon_{n+1} + \epsilon_{n+2}\\ E(y_{n+2}|y_{1:n},\theta) &=& (\phi_1^2+\phi_2)y_n + \phi_1\phi_2 y_{n-1}\\ V(y_{n+2}|y_{1:n},\theta) &=& (1+\phi_1^2)\sigma^2 \end{eqnarray*}\] Similarly, we can derive the mean and variance functions for \(y_{n+3}\) and \(y_{n+4}\): \[\begin{eqnarray*} E(y_{n+3}|y_{1:n},\theta) &=& (\phi_1^3+2\phi_1\phi_2)y_n + (\phi_1^2\phi_2 + \phi_2^2)y_{n-1}\\ V(y_{n+3}|y_{1:n},\theta) &=& (1+\phi_1^2 + (\phi_1^2+\phi_2)^2)\sigma^2\\ E(y_{n+3}|y_{1:n},\theta) &=& (\phi_1^4+3\phi_1^2\phi_2+\phi_2^2)y_n + (\phi_1^3\phi_2 + 2\phi\phi_2^2)y_{n-1}\\ V(y_{n+3}|y_{1:n},\theta) &=& (1+\phi_1^2 + (\phi_1^2+\phi_2)^2+(\phi_1^3+2\phi_1\phi_2)^2)\sigma^2 \end{eqnarray*}\]

Since all error terms, \(\epsilon_{n+h}\), for \(h=1,2,3,4\), are iid \(N(0,\sigma^2)\), the distribution of \(y_{n+h}|y_{1:n},\theta\) is also normal with mean and variance \(E(y_{n+h}|y_{1:n},\theta)\) and \(V(y_{n+h}|y_{1:n},\theta)\).

Since the above AR(2) is stationary, it follows that the limiting distribution is also Gaussian with zero mean and variance \[ Var(y_t) = \left(\frac{1-\phi_2}{1+\phi_2}\right)\{(1-\phi_2)^2-\phi_1^2)\}\sigma^2 \]

E = rep(0,4)
V = rep(0,4)
E[1] = phi1*y[n]+phi2*y[n-1]
V[1] = sig^2
E[2] = (phi1^2+phi2)*y[n]+phi1*phi2*y[n-1]
V[2] = (1+phi1^2)*sig^2
E[3] = (phi1^3+2*phi1*phi2)*y[n]+(phi1^2*phi2+phi2^2)*y[n-1]
V[3] = (1+phi1^2+(phi1^2+phi2)^2)*sig^2
E[4] = (phi1^4+3*phi1^2*phi2+phi2^2)*y[n]+
       (phi1^3*phi2+2*phi1*phi2^2)*y[n-1]
V[4] = (1+phi1^2+(phi1^2+phi2)^2+(phi1^3+2*phi1*phi2)^2)*sig^2       
       
L = min(E)-4*max(sqrt(V))
U = max(E)+4*max(sqrt(V))
xx = seq(L,U,length=1000)

c(y[n],y[n-1])

## [1] -5.487150 -3.508837

par(mfrow=c(1,1))
plot(xx,dnorm(xx,E[1],sqrt(V[1])),xlab="",ylab="Density",type="l",lwd=2)
for (h in 2:4)
  lines(xx,dnorm(xx,E[h],sqrt(V[h])),col=h,lwd=2)
points(E,rep(0,4),col=1:4,pch=16)
legend("topright",legend=c(
paste("y(",n+1,")|y(1),...,y(",n,")",sep=""),
paste("y(",n+2,")|y(1),...,y(",n,")",sep=""),
paste("y(",n+3,")|y(1),...,y(",n,")",sep=""),
paste("y(",n+4,")|y(1),...,y(",n,")",sep="")),col=1:4,pch=16,bty="n")

4 Monte Carlo-based \(h\)-steps ahead forecasting

M  = 10000   # Monte Carlo sizes
s  = 50      # Number of steps ahead
yf = matrix(0,s,M)
yf[1,] = phi1*y[n]+phi2*y[n-1]+rnorm(M,0,sig) 
yf[2,] = phi1*yf[1,]+phi2*y[n]+rnorm(M,0,sig) 
for (h in 3:s)
  yf[h,] = phi1*yf[h-1,]+phi2*yf[h-2,]+rnorm(M,0,sig) 
  
qyf = apply(yf,1,quantile,c(0.05,0.5,0.95))  
  
par(mfrow=c(2,2))
for (h in 1:4){
  xx = seq(min(yf[h,]),max(yf[h,]),length=100)
  breaks = seq(min(yf[h,]),max(yf[h,]),length=30)
  hist(yf[h,],xlab="",prob=TRUE,main="",xlim=range(xx),
       breaks=breaks)
  title(paste("y(",n+h,")|y(1),...,y(",n,")",sep=""))
  lines(xx,dnorm(xx,E[h],sqrt(V[h])),col=2,lwd=2)  
}  

c(mean(yf[s,]),sqrt(var(yf[s,])))

## [1] 0.02971917 7.12590341

sd = sqrt(sig^2*(1-phi2)/((1+phi2)*((1-phi2)^2-phi1^2)))
c(0,sd)

## [1] 0.00000 7.12069

quants = qnorm(c(0.05,0.5,0.95),0,sd)

par(mfrow=c(1,1))
range = range(y[(n-50):n],qyf)
plot((n-50):n,y[(n-50):n],lwd=3,xlim=c(n-50,n+s),ylim=range,
     xlab="Time",ylab="Observations")
points((n+1):(n+s),qyf[1,],col=3)
points((n+1):(n+s),qyf[3,],col=3)
points((n+1):(n+s),qyf[2,],col=2)
abline(h=quants[1],col=3)
abline(h=quants[2],col=2)
abline(h=quants[3],col=3)