AR(1) with missing

Data augmentation

1 AR(1) model

Let us assume that time-series observations \(y_1,\ldots,y_n\) are modeled via a autoregressive (AR) model of order one, ie. AR(1): \[ y_t|y_{t-1},\theta \sim N(\mu+\phi y_{t-1},\sigma^2), \] for \(t=2,\ldots,n\) and \(\theta=(\mu,\phi,\sigma^2)\). Inference is conditional on the value of \(y_1\) due to the lag structure of the AR(1) process. However, some of the observations are missing, as in the following simulated example, where observations \(t=65,\ldots,75\) are missing.

set.seed(31415)
mu  = 0.05
phi = 0.95
sig = 0.25
n   = 100
y   = rep(mu/(1-phi),n)
for (t in 2:n)
  y[t] = mu + phi*y[t-1]+rnorm(1,0,sig)
y.true = y
par.true = c(mu,phi,sig)

a=65
b=75
midpoint = trunc((a+b)/2)
y[a:b]=(y[a-1]+y[b+1])/2

par(mfrow=c(1,1))
ts.plot(y.true,xlab="Time",ylab="Observations",lwd=2)
lines(y,col=2,lwd=2)
abline(v=a,lty=2)
abline(v=b,lty=2)
legend("topleft",legend=c("True data","Observed data"),col=1:2,lty=1,lwd=2,bty="n")

2 Gibbs sampler

Assuming a noninformative prior for the parameters \(\mu,\phi,\sigma^2,y_a,y_{a+1},\ldots,y_b\), ie. \[ p(\mu,\phi,\sigma^2,y_a,y_{a+1},\ldots,y_b) \propto \frac{1}{\sigma^2}, \] such that posterior inference is feasible via Markov chain Monte Carlo through the Gibbs sampler. More preciselly, we can easily the full conditionals for each one of the parameters.

2.1 Full conditional of \(\mu\)

To sample \(\mu\) from its full conditional, we first define \(z_t=y_t-\phi y_{t-1}\), for \(t=2,\ldots,n\). Then, \[ p(\mu|\mbox{all}) \propto \exp\left\{-\frac{0.5}{\sigma^2} \sum_{i=1}^n (z_t-\mu)^2 \right\} \propto \exp\left\{-\frac{0.5}{\sigma^2/(n-1)} (\mu^2 - 2\mu {\bar z}) \right\}, \] where \({\bar z}=\sum_{t=2}^n (y_t-\phi y_{t-1})/(n-1)\). Therefore, the full conditional of \(\mu\) is Gaussian: \[ (\mu|\mbox{all}) \sim N({\bar z},\sigma^2/(n-1)). \]

2.2 Full conditional of \(\phi\)

To sample \(\phi\) from its full conditional, we first define \(z_t=y_t-\mu\), for \(t=2,\ldots,n\). Then, \[ p(\phi|\mbox{all}) \propto \exp\left\{-\frac{0.5}{\sigma^2} \sum_{i=1}^n (z_t-\phi y_t)^2 \right\} \propto \exp\left\{-\frac{0.5}{\sigma^2/(n-1)} \left(\phi^2\sum_{t=1}^n z_t- 2\phi\sum_{i=2}^n z_t y_{t-1}\right) \right\}. \] Therefore, the full conditional of \(\phi\) is also Gaussian: \[ (\phi|\mbox{all}) \sim N\left( \frac{\sum_{t=2}^n (y_t-\mu)y_{t-1}}{\sum_{t=2}^n y_{t-1}^2},\frac{\sigma^2}{\sum_{t=2}^n y_{t-1}^2}\right). \]

2.3 Full conditional of \(\sigma^2\)

The full conditional distribution of \(\sigma^2\) is the simplest to derive \[ p(\sigma^2|\mbox{all}) \propto \sigma^{-2} \exp\left\{-\frac{0.5}{\sigma^2} \sum_{t=2}^n (y_t-\mu-\phi y_{t-1})^2\right\}, \] which, as a function of \(\sigma^2\), is the kernal of an inverse gamma distribution with parameters \(n+1/2\) and \(\sum_{t=2}^n (y_t-\mu-\phi y_{t-1})^2/2\): \[ (\sigma^2|\mbox{all}) \sim IG\left(\frac{n+1}{2},\frac{\sum_{t=2}^n (y_t-\mu-\phi y_{t-1})^2}{2}\right) \]

2.4 Full conditional of missing \(y\)s

The full conditional for the missing \(y_t\) is given by

\[\begin{eqnarray*} p(y_t|\mbox{all}) &\propto& p(y_t|y_{t-1},\theta)p(y_{t+1}|y_t,\theta)\\ &\propto& \exp\left\{-\frac{0.5}{\sigma^2} \left[ (y_t-(\mu+\phi y_{t-1})^2+((y_{t+1}-\mu)-\phi y_t)^2 \right]\right\}\\ &\propto& \exp \left\{-\frac{0.5}{\sigma^2} \left[ y_t^2 - 2 y_t(\mu+\phi y_{t-1}) + \phi^2 y_t^2 - 2y_t\phi(y_{t+1}-\mu)\right] \right\}\\ &\propto& \exp \left\{-\frac{0.5}{\sigma^2} \left[ (1+\phi^2)y_t^2 - 2 y_t(\mu+\phi y_{t-1}+\phi(y_{t+1}-\mu))\right] \right\}, \end{eqnarray*}\] which is a Gaussian kernel with mean and variance \[ \frac{\mu(1-\phi) + \phi(y_{t-1}+y_{t+1})}{1+\phi^2} \ \ \ \mbox{and} \ \ \ \frac{\sigma^2}{(1+\phi^2)}, \] respectively.

burnin = 10000
M      = 2000
skip   = 100
niter  = burnin+M*skip
draws  = matrix(0,niter,b-a+4)
for (iter in 1:niter){
  # Sampling mu
  z       = y[2:n]-phi*y[1:(n-1)]
  mu      = rnorm(1,mean(z),sig/sqrt(n-1))
  # Sampling phi
  var.phi = sig^2/sum(y[1:(n-1)]^2)
  mean.phi = var.phi*sum((y[2:n]-mu)*y[1:(n-1)])/sig^2
  phi = rnorm(1,mean.phi,sqrt(var.phi))
  # Sampling sigma2
  par1 = n + 1
  par2 = sum((y[2:n]-mu-phi*y[1:(n-1)])^2)
  sig  = sqrt(1/rgamma(1,par1/2,par2/2))
  # Sampling missing observations
  for (t in a:b)
    y[t] = rnorm(1,(mu*(1-phi)+phi*(y[t-1]+y[t+1]))/(1+phi^2),sig/sqrt(1+phi^2))
  draws[iter,] = c(y[a:b],phi,sig,mu)
}
ind  = seq(burnin+1,niter,by=skip)
mus  = draws[ind,b-a+4]
phis = draws[ind,b-a+2]
sigs = draws[ind,b-a+3]
ymid = draws[ind,midpoint-a+1]

2.5 Posterior summaries

par(mfrow=c(3,4))
ts.plot(mus,xlab="Iterations",ylab="",main=expression(mu))
ts.plot(phis,xlab="Iterations",ylab="",main=expression(phi))
ts.plot(sigs,xlab="Iterations",ylab="",main=expression(sigma))
ts.plot(ymid,xlab="Iterations",ylab="",main=paste("y[",midpoint,"]",sep=""))
acf(mus,main="")
acf(phis,main="")
acf(sigs,main="")
acf(ymid,main="")
hist(mus,prob=TRUE,xlab="",main="");abline(v=par.true[1],col=2,lwd=2)
hist(phis,prob=TRUE,xlab="",main="");abline(v=par.true[2],col=2,lwd=2)
hist(sigs,prob=TRUE,xlab="",main="");abline(v=par.true[3],col=2,lwd=2)
hist(ymid,prob=TRUE,xlab="",main="");abline(v=y.true[midpoint],col=2,lwd=2)

2.6 Posterior for the missing points

ys = draws[,1:(b-a+1)]
qy = t(apply(ys,2,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(1,1))
plot(y.true,ylim=range(y.true,draws[,1:(b-a+1)]),xlab="Time",ylab="Observations")
abline(v=a,lty=2)
abline(v=b,lty=2)
for (i in 1:M)
  lines(a:b,draws[i,1:(b-a+1)],col=grey(0.8))
points(a:b,y.true[a:b],pch=16)
lines(a:b,qy[,1],col=2,lwd=2)
points(a:b,qy[,2],col=2,pch=16)
lines(a:b,qy[,3],col=2,lwd=2)
lines(y.true)