2
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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... . $\endgroup$
    – whuber
    Commented Dec 14, 2018 at 16:33
  • $\begingroup$ Fair enough, the proper assumption would be Lognormal I guess? I will adapt the question, thanks :-) $\endgroup$
    – germannp
    Commented Dec 16, 2018 at 11:04
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Dec 16, 2018 at 20:51
  • $\begingroup$ 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. $\endgroup$
    – germannp
    Commented Dec 17, 2018 at 8:54
  • $\begingroup$ 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. $\endgroup$
    – germannp
    Commented Dec 20, 2018 at 10:39

2 Answers 2

1
$\begingroup$

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
}
$\endgroup$
2
  • $\begingroup$ 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 ... $\endgroup$
    – germannp
    Commented Dec 12, 2018 at 8:57
  • $\begingroup$ 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. $\endgroup$ Commented Dec 12, 2018 at 19:54
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.