# Timeseries with multiplicative noise in Stan

Say we have a monthly time series $$y_t \geq 0$$ dominated by seasonality, where the absolute differences from year to year are much smaller during low season. To avoid negative values and capture the higher uncertainty during high seasons I want to fit the time series based on last year's value with multiplicative noise, i.e. as

$$y_t = y_{t-12} \cdot \varepsilon_t.$$

What is the correct Stan code for this?

• What is your argument that $y_t$ must have a Normal distribution? After all, if that is correct, then it applies just as well to $y_{t-12},$ implying $y_t$ is the product of two independent Normally distributed variables, which definitely is not Normal. Something's not quite right here... .
– whuber
Commented Dec 14, 2018 at 16:33
• Fair enough, the proper assumption would be Lognormal I guess? I will adapt the question, thanks :-) Commented Dec 16, 2018 at 11:04
• Could you explain what the data are supposed to be? You seem be generating them as a perturbation of a sine wave, but that would have a standard deviation different from that of the perturbing lognormal noise. Moreover, although $\sigma$ might equal $0.1,$ that is not the same as the standard deviation of that lognormal variable--it would be the geometric sd.
– whuber
Commented Dec 16, 2018 at 20:51
• The data is supposed to simulate freight train traffic volumes. They show some trends, but are mainly periodic. Some vary quite a bit during high seasons, e.g. crops are not always harvested the same week -- but for the rest there is none of this crop and therefore no traffic with high certainty. Commented Dec 17, 2018 at 8:54
• I made a mistake simulating the data in the current version of the question ... I will edit the question heavily and adapt the current version as a reply, to make it useful for future reference. Commented Dec 20, 2018 at 10:39

I would try exponential noise, in which case your Stan program would be

data {
int<lower = 0> N;
vector<lower = 0>[N] y;
}
transformed data {
vector[N - 12] y_lead = y[13:N];
vector[N - 12] y_lag_inv = inv(y[1:(N - 12)]);
}
parameters {
real<lower = 0> alpha;
}
model {
// prior on alpha
}

• Thanks for your input. Would you mind explaining the rational behind exponential noise? I believe exponentially distributed numbers are positive and mainly around zero, so I would have to add them when generating predictions? I clarified the question. Also, what does inv() do? I can't find it in the manual ... Commented Dec 12, 2018 at 8:57
• inv() is the inverse, i.e. reciprocal. The expectation of this exponential distribution is alpha * y_lag so it is a reasonable starting point. It is skewed toward zero, which may be unrealistic. If so, you could then move to a gamma distribution and estimate an additional parameter. Commented Dec 12, 2018 at 19:54

One might argue, that multiplicative errors should follow the log-normal distribution, because a log-normally distributed value times another one is again log-normally distributed. As the median of $$\textrm{Lognormal}(\mu, \sigma)$$ is $$e^\mu$$, the model

$$y_t \sim \textrm{Lognormal}(\log(y_{t-12}), \,\sigma)$$

might be reasonable. This can be modeled in Pystan as:

import numpy as np
import pystan

n_years = 4
y = (np.sin(np.linspace(0, n_years*np.pi, n_years*12+1))**2  + 0.05) * np.random.lognormal(0, 0.1, n_years*12+1)

model = pystan.StanModel(model_code="""
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real<lower=0, upper=0.5> sigma;
}
model {
y[13:N] ~ lognormal(log(y[1:N-12]), sigma);
}
""")

fit = model.sampling(data={'N': len(y), 'y': y}, iter=1000, chains=4)


Which recovers $$\sigma$$ fairly well:

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
sigma   0.15  7.1e-4   0.02   0.12   0.14   0.15   0.16   0.19    704    1.0
lp__   49.62    0.03   0.76  47.41  49.46   49.9   50.1  50.17    860   1.01


The fit can be further improved by using the (geometric) mean of the last two year's values.