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.).

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

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).

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:N1], sigma_mu);
for(i in 365:N)
period_yearly[i] ~ normal(sum(period_yearly[i364:i1]), 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.