Time series model

Interfaces RStan
1

Hi there,
I’m new to coding time series models with Stan. I’ve got below R code to generate the data and a Stan model that fits the data. I think it’s all ok, but it would be great if someone who’s worked with time series models before could have a look.
I’m not quite sure whether i’ve done the priors correctly.

R code to generate the data:
Each outcome is drawn from a normal distribution of mean my.mus for that trial (with sd=my.sd).
These means also drift over time, following a normal distribution (with noise my.k).

my.k=19
my.sd=1
my.mu1=50
ntr=20
my.outcomes = array(NA,dim=ntr)
my.mus      = array(NA,dim=ntr)
my.mus[1]   = my.mu1
for (itr in 1:ntr){
  my.outcomes[itr] = round(rnorm(1,my.mus[itr],my.sd));
  if (itr<ntr){
    my.mus[itr+1] = rnorm(1,my.mus[itr],my.k);
  }
}

From looking at the Stan manual, this is not quite an autoregressive model I think as the outcome on the next trial does not depend on the outcome of the last trial, but instead there is a relationship over time between the (unobserved) means

Here is the Stan code:

data {
  int<lower=0> ntr;
  int outcomes[ntr];
}

parameters {
  real mag_mean[ntr];
  real<lower=0> mag_sd;
  real<lower=0> k;
}

model {
  // Priors
  mag_mean[1] ~ normal(50,20); // priors for other trials are given in loop (?)
  mag_sd      ~ normal(0,5);
  k           ~ normal(0,5);
  
  // Model
  for (itr in 1:(ntr)){
    outcomes[itr]     ~ normal(mag_mean[itr],mag_sd);
    if (itr<ntr){
      mag_mean[itr+1]   ~ normal(mag_mean[itr],k);
    }
  }
}

Many thanks
Jacquie

2

Are there multiple observations at each time point? If not, then your priors for the variance parameters are going to be doing a lot of work to render the model identifiable. Have you considered a Gaussian Process model?

3

Actually there was a previous post on a similar topic, sorry I had missed that before:

which seems to tell me that I should vectorise and use the ‘matt trick’ … trying this now

4

Thanks for the quick reply.
I’m looking at the Gaussian Process model chapter section 10.1 in current manual. It seems that these models have y and x? But my model only has observations y, nothing else is observed. Am I misunderstanding these models? (Well possible as I did not understand the text well at all given my lack of knowledge about this)

5

Why dont you write it as Kalman Filter? This Handook could help you alot. And with the gaussian_dlm_obs_lpdf function is really ease to write your model. I think ctsem package could do it for you. But the right person to ask is @Charles_Driver.

6

In the SUG, x is analogous to time. There’s a worked-through example of a GP regression (with inference on the latent means, which I suspect is what you’re after) here.

You should also take a look at the lgpr for a possibly more user-friendly interface to GPs.

1 Like