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?
 A: 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 {
  y_lead ~ exponential(inv(alpha) * y_lag_inv);
  // prior on alpha
}

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