Structural time series with seasonality, sampling very slowly

[scellus]

The whitened model below mixes well with default sampler parameters, 1000 iterations in 168 seconds.

(This was an interesting experiment in that I have my own AR model with gaussian likelihood but lots of missing data. It would benefit of the same treatment. ;))

data {
  int<lower=1> n;   // n = 2632
  int y[n]; // number of shootings per day
}
transformed data {
  int<lower=1,upper=365> yday[n];
  for (i in 1:n) yday[i] = i % 365 + 1;
}
parameters {
  vector[n] mu_innovations;
  vector[365] seasonal_innovations; // yearly seasonality 
  real<lower=0> sigma_mu;
  real<lower=0> sigma_yday;
  real baseline;
}
transformed parameters {
  // zero-mean seasonal term
  vector[365] y_seasonal;
  vector[n] mu;
  { vector[365] seasonal_with_trend;
    real trend;
    seasonal_with_trend = cumulative_sum(seasonal_innovations);
    trend = seasonal_with_trend[365];
    for (i in 1:365)
        y_seasonal[i] = sigma_yday/100 * (seasonal_with_trend[i] - trend * i/365.0); }

  mu = sigma_mu/100 * cumulative_sum(mu_innovations);
}
model {
  seasonal_innovations ~ normal(0, 1);
  mu_innovations ~ normal(0, 1);
  y ~ poisson_log(baseline + mu + y_seasonal[yday]);
  sigma_mu ~ lognormal(-3.5 + log(100), 2);
  sigma_yday ~ lognormal(-3.5 + log(100), 2);
}