Time series random walk parameters

If I have a time series generated by a additive random walk

y_n = y_{n-1} - R + \sigma

If I know R can only be positive, how can I model this to get the best possible estimate for R, which in this case is a real quantity of interest. This is not a predictive exercise, I want to know what R was.

I’ve found with a model such as:

parameters {
  real<lower=0> sigma;
  real<lower=0> R;
model {
  R ~ normal(0, 1;
  sigma ~ cauchy(0, 1);
  y[2:N] ~ normal(y[1:(N-1)] - R, sigma)

That R and sigma tend to compete, sigma is overestimated, R is underestimated.

Is this just not something I can model well?

i’m generated my data by:

R = 5
y = vector(length=500)
y[1] = 100
sigma = 1
for(i in 2:length(y)){
    y[i] = y[i-1] - R + norm(1, 0, sigma)

i should say that this is the simplest version of the problem I am looking at.
The true time series would be something like

y_n = y_{n-1} - R + M + J + \sigma

where M and J have uncertainties and time dependence.


1 Like

They’re fine alternative explanations for the pattern because the normal distribution is agnostic w.r.t. the balance of l/r moves so they will compete. A model that might work better is to regress y on t with R as the slope.

You can add a parameter on y_{n-1} that will act as a weight for y_{n-1} versus the other covariates, i.e.

y_{n} \sim N(\beta y_{n-1} + E,\sigma)

\beta should be constrained to the (-1,1) interval, either with a prior around 0 if you have plenty of data or with a hard constraint if you don’t. Otherwise you could end up with an explosive time series estimate.

Thanks, so like an AR(1)?

I’ve just tried it with y_n \sim N(\beta y_{n-1} - R, \sigma).
First with hard bounds [-1, 1], which produces divergent transitions, and with all of the mass concentrated on 1.
with a no bounds and \beta \sim N(0, 1) I get a distribution of \beta that has most of it’s mass > 1 (mean= 1.0015). both of these produce worse estimates for R.

Have i missed something from your suggestion?

Sorry about that, I forgot that you have actual data (I usually use these in latent variable models).

Probably all you need to do is switch to a non-centered parameterization (similar to how you generated the data). That often does dramatically better with sampling in time series models. Here’s some sample code:

If that doesn’t work and you want to switch to a stationary model, then the most straightforward thing is to difference your data to make it stationary:

\Delta y_n = \beta \Delta y_{n-1} - E + \epsilon_n

This is a slightly different model, especially if you don’t difference E. But generally it’s easier to work with stationary time series models.

I still think I must be missing something here, sorry.

But won’t transforming my data to be stationary remove the signal I am interested in? It’s the non-stationary R I’m trying to estimate. (I presume your E = my R)

Similarly, I have no sampling problems as far as I can tell, the model converges quickly, is there another benefit of a non-centred param I’m unaware of?

You’ll need to update the toy DGP to match your model - the true value is on the parameter boundary hence the divergent transitions. e.g. try this with your AR(1) code instead:
y[i] = 0.7 * y[i-1] - R + rnorm(1, 0, sigma)

Not sure what your real data looks like though - is a random walk i.e. explosive AR(1) model reasonable/within physical limitations? If so, scrap the autoregressive parameter and the first differences of y give you a constant mean model (about your parameter of interest), which will be extremely fast to estimate:

data {
  int<lower=1> N;
  vector[N] y;
transformed data {
  vector[N-1] ydiff = y[2:N] - y[1:(N-1)];
parameters {
  real<lower=0> sigma;
  real<lower=0> R;
model {
  R ~ normal(0, 1);
  sigma ~ cauchy(0, 1);
  ydiff ~ normal( -R, sigma);

Not sure of what is intended with the M/J terms later but you might want to explore a state-space representation as complexity grows - see e.g. here https://github.com/jrnold/ssmodels-in-stan