I wish to include holidays/events(similar to FBprophet) to stan model

I wrote a code using MCMC algorithm where it could calculate the seasonality of a year(Please correct me if my logic and understanding about the algorithm is wrong since I am new to it.).

  1. I wish to add the same time extract the weekly covariance and also the weekly(vectors).

  2. Is there a possible way to implement it the same as FBProphet where it could also input the holidays and events to the stan code? - I tried to check from many sites but could not find an answer and I finally arrived to this site and hope that all the experts here could assist me in editing my code so that it could include the event/holiday(possible ways to extract the information as FBProhphet does).-

  3. If possible, could you provide some possible reference and also some better way of including more regressors to the stan code? In FBProphet we could easily use add_regressor.

Even if just one of the above could be shared is greatly appreciated too.

Below is my stan code

data {
    int N;
    int N_pred;
    vector[N] Y;

parameters {
    real<lower=0> sigma_mu;
    real<lower=0> sigma_period_yearly;
    real<lower=0> sigma_Y;
    real c_ar[2];
	real mu0;
    vector[N] mu;
    vector[N] period_yearly;


transformed parameters {
    vector[N] y_mean;
    y_mean = mu + period_yearly;

model {
    mu[2:N] ~ normal(mu[1:N-1], sigma_mu);
    for(i in 365:N)
      period_yearly[i] ~ normal(-sum(period_yearly[i-364:i-1]), sigma_period_yearly);
    Y ~ normal(y_mean, sigma_Y);

generated quantities {
    vector[N + N_pred] mu_all;
    vector[N + N_pred] period_yearly_all;
    vector[N_pred] y_new;
    mu_all[1:N] = mu;
    period_yearly_all[1:N] = period_yearly;
    for (t in 1:N_pred) {
        mu_all[N + t] = normal_rng(mu_all[N + t - 1], sigma_mu);
        period_yearly_all[N + t] = normal_rng(-sum(period_yearly_all[(N + t - 11):(N + t - 1)]), sigma_period_yearly);
        y_new[t] = normal_rng((mu_all[N + t] + period_yearly_all[N + t]	) + trend_all[N + t] + ar_all[N + t], sigma_Y);

Thank you and please kindly assist me.

Is there any experts who could assist me in this please? If not, just an example of ideas of how to implement it. I have tried to search and learn from many other places for 3 days now and still could not understand how to implement it.

You help is kindly appreciated.

Hi! Well, an expert will be @Charles_Driver of cstem package, I think he have deal with models like yours before. And you should check this post it talks of something similar to yours.

Another expert will be Sean Taylor of prophet but I don’t know if he is here, you can see prophet’s Stan code here. But prophet use Generalized additive Models so I think it might be a different model.

For multiple periodicities you can use BRMS and Fourier bases. If I am not wrong @avehtari and @paul.buerkner do it with GP in their article. But check it maybe helps. The used code is here.

Another approach of before (I think) will be using Fourier decomposition base, fourier() function from forecast package does it for you, It can handle multiple periodicities too and you can fit your model in Stan like a simple regression. An idea of it is in Rob’s online book

I hope this works for now! I will try to check your code today. :)


@asael_am Thank you so much for sharing some of the topics and codes that had been discussed. I read through a few of the articles and try to understand and compare the codes and you are right your first post is related to what I try to deal with(seasonality by week but they did not specify how they conduct it. I will try to ask them to share their findings.)

I also attempted to read from other articles in which you provided but I still do not know how to include the holidays/events as per FBProphet could as well as including additional regressors. Thank you again for trying to review my code and I will be looking forward to your suggestions.