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