Why are my predictied values from a Bayesian AR(1) model lagging behind the data?

Question

Summary: I have simulated some data on an AR(1) process in R and fit the model in Stan. When plotting the predictions, the predicted values tend to lag behind the true values. Why is this?

Detail

I am simulating fake data from the following AR(1) model:

\begin{equation} y_{t} \sim Normal(\eta_{t}, \sigma) \end{equation}

\begin{equation} \eta_{t} = \mu + \alpha y_{t-1} \end{equation}

for $t$ from $2,...,N$, and where $y_{1} \sim Normal(\mu, 0.1)$

In R:

# requirements
require(tidyverse)
require(coda)
require(rstan)

# convenience function for the highest density interval
hdi <- function(x, bound=NULL){
  hdi_ <- coda::HPDinterval(coda::as.mcmc(x))
  
  if(bound=="low"){
    hdi_ <- hdi_[[1]]
  }
  if(bound=="high"){
    hdi_ <- hdi_[[2]]
  }
  
  return(hdi_)
}

set.seed(2021)

# number of data points
N <- 30
# mean
mu <- 0.5
# AR(1) coefficient
alpha <- 0.3
# noise
sigma <- 0.1

# simulate y
y <- numeric(N)
y[1] <- rnorm(1, mu, 0.1)
for(t in 2:N){
  y[t] = rnorm(1, mu + alpha*y[t-1], sigma)
}

# package data
d <- data.frame(t=1:N, y=y)

I then fit the model to this data using the following Stan program to return the parameter values:

data{
  int<lower=0> N;
  vector[N] y;
}

parameters{
  real mu;
  real alpha;
  real<lower=0> sigma;
}

transformed parameters{
  vector[N-1] eta;
  eta = mu + alpha * y[1:(N-1)];
}

model{
  mu ~ normal(0, 1);
  sigma ~ normal(0, 1);
  alpha ~ normal(0, 1);
  
  y[2:N] ~ normal(eta, sigma);
}

generated quantities{
  vector[N-1] log_lik;
  vector[N-1] PPD;
  for(t in 1:(N-1)){
    log_lik[t] = normal_lpdf(y[t] | eta[t], sigma);
    PPD[t] = normal_rng(eta[t], sigma);
  }
}

Fit the model:

# compile the Stan model
ar1 <- rstan::stan_model("ar1.stan")

# fit the model
fit_ar1 <- rstan::sampling(
  ar1, 
  data = list(N=N, y=d$y), 
  cores = 4
)

# posterior
ps <- as.data.frame(fit_ar1)

# package the data and predictions into a df
pred_df <- data.frame(
  t = 1:N,
  y = d$y, 
  pred_mean = c(NA, (apply(ps[,grep("eta", colnames(ps))],2,mean))),
  pred_hdi_low = c(
    NA, 
    (apply(ps[,grep("eta", colnames(ps))],2, function(x) hdi(x, "low")))
  ),
  pred_hdi_high = c(
    NA, 
    (apply(ps[,grep("eta", colnames(ps))],2, function(x) hdi(x, "high"))) 
  ),
  ppd_hdi_low = c(
    NA, 
    (apply(ps[,grep("PPD", colnames(ps))],2, function(x) hdi(x, "low"))) 
  ),
  ppd_hdi_high = c(
    NA, 
    (apply(ps[,grep("PPD", colnames(ps))],2, function(x) hdi(x, "high"))) 
  )
)

# plot the predictions against the data
ggplot(pred_df) +
  geom_ribbon(aes(x=t, ymin=pred_hdi_low, ymax=pred_hdi_high), 
              fill=alpha("green", 0.7)) + 
  geom_ribbon(aes(x=t, ymin=ppd_hdi_low, ymax=ppd_hdi_high), 
              fill=alpha("yellow", 0.4)) + 
  geom_line(aes(t, y)) + 
  geom_point(aes(t, y)) + 
  geom_line(aes(t, pred_mean), colour="darkolivegreen", size=1) + 
  geom_point(aes(t, pred_mean), colour="darkolivegreen") + 
  theme(panel.grid = element_blank()) + 
  scale_x_continuous(limits=c(2, N))

In this plot, the black line and points show the raw data, the dark green line and ribbon the posterior mean and uncertainty around the mean, and the yellow ribbon the 95% HDI of the posterior predictive distribution. But it seems the predicted values tend to lag behind the real values, in this case by one data point. An increase in the real data is reflected by a increase in the predicted value shifted by one time point in the future.

This is a recurring theme for my exploration of time series models, even on data that have been simulated from that model as I have done here.

Does anyone have any advice?

Answer 1

I think what you see is mostly caused by the forecast $\hat y_t$ being a function of $y_{t-1}$ only which (sometimes but not always, see below) makes $\hat y_t$ more similar to $y_{t-1}$ than $y_t$.

If we for simplicity ignore the uncertainty in the parameters, the one step ahead predicted values based on all data up to time $t-1$ for an AR(1) model $y_t=\phi y_{t-1}+w_t$ are given by $$ \hat y_t=E(y_t|y_{t-1},y_{t-2},\dots)=\phi y_{t-1}. $$ The mean square errors quantifying the similarity of $\hat y_t$ and, respectively, $y_t$ and $y_{t-1}$, are $$ E((\hat y_t-y_t)^2)=\mbox{Var}(\phi y_{t-1} - y_t)=\mbox{Var}(-w_t)=\sigma^2 $$ and \begin{align} E((\hat y_t-y_{t-1})^2)&=\mbox{Var}(\phi y_{t-1} - y_{t-1})\\&=(1-\phi)^2\mbox{Var}(y_{t-1}) \\&=(1-\phi)^2\frac{\sigma^2}{1-\phi^2}\\&=\sigma^2\frac{1-\phi}{1+\phi}. \end{align}

The first mean square error is larger than the second one (the effect you see in your data) when $$ 1>\frac{1-\phi}{1+\phi} $$ or $$ 1+\phi>1-\phi $$ or $$ \phi>0. $$