1
$\begingroup$

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))

Predictions

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?

New contributor
user_15 is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
0
$\begingroup$

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. $$

$\endgroup$
  • $\begingroup$ Thanks, I'm starting to understand. What I still don't understand is that the posterior predictive distribution should match the data, i.e. simulating new observations should match the original data, but they are also out of sync. $\endgroup$ – user_15 Jan 14 at 10:09
  • $\begingroup$ I realised that my posterior predictions were not actually taking into account the stochastic process correctly. In producing PPD, I should have used previous values of PPD to construct a new time series. However, I just used the values of eta directly. $\endgroup$ – user_15 Jan 14 at 16:56

Your Answer

user_15 is a new contributor. Be nice, and check out our Code of Conduct.

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.