Holt time-series forecasting in Stan

holt.stan

data {
  int<lower=3> n;
  vector[n] y;
  int<lower=0> h;
}
parameters {
  real<lower=0, upper=1> alpha;
  real<lower=0, upper=1> beta;
  real<lower=0> sigma;
  vector[2] mu_init;
}
transformed parameters {
  real l;
  real lp;
  real b;
  vector[n] mu;

  lp = y[1];
  l = lp;
  b = y[2] - y[1];
  mu[1:2] = mu_init;
  
  for (t in 2:(n-1)) {
    l = alpha * y[t] + (1 - alpha) * (lp + b);
    b = beta * (l - lp) + (1 - beta) * b;
    mu[t+1] = l + b;
    lp = l;
  }
}
model {
  y ~ normal(mu, sigma);
}
generated quantities {
  vector[h] yhat;
  
  for (k in 1:h)
    yhat[k] = normal_rng(l + b*k, sigma);
}