Forecast method that allows a continuous assessment of the forecast likelihood over a horizon?

Hello there, I’m very new to Bayesian, let alone Stan. I recently came across a piece in a forecast report that looks like this:

forecast

In essence, this forecast allows an estimate of forecast uncertainty represented by the forecast intervals. More importantly, it allows a continuous assessment of the likelihood of the outcome at the end of the forecast horizon based on the update of the data. For example, it can deliver messages such as: Conditional on the data updated at Oct. 2022, the probability that the sum of actual values is greater than that of forecast values by the end of the forecast horizon is 64%.

The document doesn’t describe the method, so I wonder if Bayesian modelling would be a potential solution. My first question is: Conceptually, is Bayesian a decent approach to enable such a continuous assessment of probability distribution of the forecast? If yes, my plan is:

(1) Generate distributions over the forecast horizon
(2) Use the mean values of these distributions to form the forecast series
(3) Generate the intervals based on the forecast series

But I’m still not clear how to construct the distribution. The new data will allow an update of the posterior distributions, but I’m stuck going from there to the distribution of the sum over the forecast horizon.

Below is a simple Bayesian model I’m practicing with. At the moment I’m not concerned too much with details of model specification, but instead seeking advice on the conceptual direction. It’d be great if you could give me a few keywords, concepts, papers or book chapters to look into.

stan_code_simple = """
data {
  int<lower=1> T;
  int<lower=1> T_forecast;
  vector[T] y; 
}

parameters {
  real<lower=0> sigma_mu; 
  real<lower=0> sigma_y;
  vector[T] mu_err;
}

transformed parameters {
  vector[T] mu; // local trend
  mu[1] = mu_err[1];
  for(t in 2:T){
    mu[t] = mu[t-1] + mu_err[t];
  }  
}

model { 
// Priors
  sigma_mu ~ uniform(0, 20);
  sigma_y ~ uniform(0,20);
  mu_err ~ normal(0,sigma_mu);
  for(t in 1:T){
    y[t] ~ normal(mu[t], sigma_y);
  }
}

// Forecasts 
generated quantities {
  real y_forecast[T_forecast];
  real mu_forecast[T_forecast];
  mu_forecast[1] = normal_rng(mu[T], sigma_mu);
  for (t in 2:T_forecast) {
    mu_forecast[t] = normal_rng(mu_forecast[t-1], sigma_mu);
  }
  for (t in 1:T_forecast) {
    y_forecast[t] = normal_rng(mu_forecast[t], sigma_y);
  }
}
"""

Note: This question is partially cross-posted on stackexchange.com