De l’approche bayésienne au modèle de série temporelle

Modèle BSTS pour la prédiction du cours de l’or

Introduction

Le modèle BSTS (Bayesian Structural Time Series) est une approche puissante pour la modélisation et la prédiction de séries temporelles, combinant des composantes structurelles (tendance, saisonnalité, régresseurs) avec une inférence bayésienne. Ce document détaille son application à la prédiction du cours de l’or, en explicitant la formulation mathématique, l’inférence bayésienne et le choix des distributions a priori.

---

1. Formulation mathématique du modèle BSTS

Le modèle BSTS décompose la série temporelley_t(ici, le cours de l’or) en plusieurs composantes additives :

y_{t} = \underbrace{\mu_{t}}_{trend} + \underbrace{\gamma_{t}}_{seasonal} + \underbrace{\beta^{T}x_{t}}_{regression} + \epsilon_{t}

• y_{t}: cours de l’or au temps t.
  
• \mu_{t}: tendance (ex : hausse constante).
  
• \gamma_{t}: saisonnalité (ex : effet du Nouvel An chinois).
  
• \beta^{T}x_{t}: effets des variables explicatives (ex : taux directeur de la FED).
  
• \epsilon_{t}: bruit blanc gaussien,\epsilon_t \sim \mathcal{N}(0, \sigma^2_\epsilon).

---

1.1. Modélisation de la tendance locale linéaire

La tendance est modélisée par un modèle local linéaire :

\begin{aligned} \mu_{t} &= \mu_{t-1} + \delta_{t-1} + u_{t}, \quad u_t \sim \mathcal{N}(0, \sigma^2_u) \\ \delta_{t} &= \delta_{t-1} + v_{t}, \quad v_t \sim \mathcal{N}(0, \sigma^2_v) \end{aligned}

• \mu_t: niveau de la série au temps t.
  
• \delta_t: pente de la tendance.
  
• u_t, v_t: bruits gaussiens de variances\sigma^2_u, \sigma^2_v.

---

1.2. Modélisation de la saisonnalité stochastique

La saisonnalité est modélisée par une somme de sinusoïdes stochastiques :

\gamma_{t} = -\sum_{s=1}^{S-1}\gamma_{t-s} + w_{t}, \quad w_t \sim \mathcal{N}(0, \sigma^2_w)

• S: période saisonnière (ex : 12 pour des données mensuelles).
  
• w_t: bruit saisonnier, de variance\sigma^2_w.

---

2. Inférence bayésienne dans BSTS

L’inférence bayésienne consiste à déduire la distribution postérieure des paramètres p(\theta|X,y), où \theta représente l’ensemble des paramètres du modèle,X la matrice des régresseurs etyle vecteur des observations.

2.1. Théorème de Bayes

\underbrace{p(\theta|X,y)}_{postérieur} \propto \underbrace{p(y|X,\theta)}_{vraissemblance} \times \underbrace{p(\theta)}_{prior}

• Vraissemblance :p(y|X,\theta)mesure la compatibilité des données avec les paramètres.
  
• Prior :p(\theta)encode les connaissances a priori sur les paramètres.

2.2 Exemple illustratif en R

library(statsr)
library(ggplot2)
library(tidyverse)

data(tapwater)
glimpse(tapwater)

Rows: 28
Columns: 6
$ date       <fct> 2009-02-25, 2008-12-22, 2008-09-25, 2008-05-14, 2008-04-14,…
$ tthm       <dbl> 34.38, 39.33, 108.63, 88.00, 81.00, 49.25, 75.00, 82.86, 85…
$ samples    <int> 8, 9, 8, 8, 2, 8, 6, 7, 8, 4, 4, 4, 4, 6, 4, 8, 10, 10, 10,…
$ nondetects <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ min        <dbl> 32.00, 31.00, 85.00, 75.00, 81.00, 26.00, 70.00, 70.00, 80.…
$ max        <dbl> 39.00, 46.00, 120.00, 94.00, 81.00, 68.00, 80.00, 90.00, 90…

ggplot(tapwater, aes(x = tthm)) +
  geom_histogram(bins = 10, fill = "#FF6B6B", alpha = 0.7) +
  labs(
    title = "Histogramme de pph",
    x = "pph",
    y = "Fréquence"
  ) +
  theme_minimal()

set.seed(42)

# Hyperparamètres du prior
m_0 = 35
n_0 = 25
s2_0 = 156.25
v_0 = n_0 - 1

# Statistiques de l'échantillon
Y = tapwater$tthm
ybar = mean(Y)
s2 = var(Y)
n = length(Y)

# Hyperparamètres du posterior
n_n = n_0 + n
m_n = (n*ybar + n_0*m_0) / n_n
v_n = v_0 + n
s2_n = ((n-1)*s2 + v_0*s2_0 + n_0*n*(m_0 - ybar)^2/n_n) / v_n

# Inference statsr
bayes_inference(
  tthm, data = tapwater, prior = "NG",
  mu_0 = m_0, n_0 = n_0, s_0 = sqrt(s2_0), v_0 = v_0,
  stat = "mean", type = "ci", method = "theoretical",
  show_res = TRUE, show_summ = TRUE, show_plot = FALSE
)

Single numerical variable
n = 28, y-bar = 55.5239, s = 23.254
(Assuming proper prior:  mu | sigma^2 ~ N(35, *sigma^2/25)
(Assuming proper prior: 1/sigma^2 ~ G(24/2,156.25*24/2)

Joint Posterior Distribution for mu and 1/sigma^2:
 N(45.8428, sigma^2/53) G(52/2, 8.6769*52/2)

Marginal Posterior for mu:
Student t with posterior mean = 45.8428, posterior scale = 2.9457 on 52 df

95% CI: (39.9319 , 51.7537)

# -----------------------------
# 1. Distribution de sigma (prior & posterior)
# -----------------------------

# Prior : phi = 1/sigma^2 ~ Gamma(v0/2, v0*s2_0/2)
phi_prior = rgamma(5000, shape = v_0/2, rate = v_0*s2_0/2)
sigma_prior = 1/sqrt(phi_prior)

# Posterior : phi_n = 1/sigma^2 ~ Gamma(vn/2, vn*s2_n/2)
phi_post = rgamma(5000, shape = v_n/2, rate = v_n*s2_n/2)
sigma_post = 1/sqrt(phi_post)

df_sigma = tibble(
  sigma = c(sigma_prior, sigma_post),
  type = rep(c("Prior", "Posterior"), each = 5000)
)

ggplot(df_sigma, aes(x = sigma, fill = type, color = type)) +
  geom_density(alpha = 0.3, linewidth = 1) +
  labs(
    title = "Distribution a priori et a posteriori de σ",
    x = "σ",
    y = "Densité"
  ) +
  theme_minimal()

# -----------------------------
# 2. Distribution de mu (prior & posterior)
# -----------------------------

# Prior : mu ~ Normal(m0, sigma/sqrt(n0))
mu_prior = rnorm(5000, mean = m_0, sd = sigma_prior / sqrt(n_0))

# Posterior : mu_n ~ Normal(mn, sigma/sqrt(nn))
mu_post = rnorm(5000, mean = m_n, sd = sigma_post / sqrt(n_n))

df_mu = tibble(
  mu = c(mu_prior, mu_post),
  type = rep(c("Prior", "Posterior"), each = 5000)
)

ggplot(df_mu, aes(x = mu, fill = type, color = type)) +
  geom_density(alpha = 0.3, linewidth = 1) +
  labs(
    title = "Distribution a priori et a posteriori de μ",
    x = "μ",
    y = "Densité"
  ) +
  theme_minimal()

---

2.2. Choix des distributions a priori

Chaque composante du modèle BSTS a ses propres paramètres, auxquels on assigne des priors Inverse-Gamma pour les variances et des priors Normaux pour les coefficients.

a. Tendance

Paramètre	Prior typique	Rôle
\sigma^2_u	Inv-Gamma(0.01, 0.01)	Variabilité du niveau de la tendance.
\sigma^2_v	Inv-Gamma(0.01, 0.01)	Variabilité de la pente de la tendance.

b. Saisonnalité

Paramètre	Prior typique	Rôle
\sigma^2_w	Inv-Gamma(0.1, 0.1)	Variabilité des amplitudes saisonnières.

c. Régresseurs

Paramètre	Prior typique	Rôle
\beta_j	Normal(0,\sigma^2_\beta)	Coefficients de régression.
\sigma^2_\beta	Inv-Gamma(0.01, 0.01)	Variabilité des coefficients.

d. Bruit global

Paramètre	Prior typique	Rôle
\sigma^2_\epsilon	Inv-Gamma(0.001, 0.001)	Variabilité résiduelle non expliquée.

---

3. Implémentation pratique avec bsts en R

Voici un exemple minimal d’implémentation du modèle BSTS pour prédire le cours de l’or :

library(bsts)
library(readr)
library(dplyr)
library(lubridate)
library(tidyverse)

path <- getwd()
gold <- read.csv(paste0(path,"/data/monthly.csv"))

gold <- gold %>%
  mutate(
    Date = paste0(Date, "-01"),  # transforme 2025-04 → 2025-04-01
    Date = as.Date(Date)
  ) %>%
  arrange(Date)

ggplot(gold, aes(x = Date, y = Price)) +
  geom_line(color = "#0072B2", linewidth = 1) +
  labs(
    title = "Cours de l'or en dollar",
    x = "Mois",
    y = "Prix"
  ) +
  theme_minimal()

#Modelisation
y <- gold$Price %>% as.double()
# Mise en log
y <- log(y)

ss <- AddLocalLinearTrend(list(), y)    
ss <- AddSeasonal(ss, y, nseasons = 12)
model1 <- bsts(y,state.specification = ss,niter = 1000)

=-=-=-=-= Iteration 0 Fri Apr  3 14:12:22 2026 =-=-=-=-=
=-=-=-=-= Iteration 100 Fri Apr  3 14:12:25 2026 =-=-=-=-=
=-=-=-=-= Iteration 200 Fri Apr  3 14:12:27 2026 =-=-=-=-=
=-=-=-=-= Iteration 300 Fri Apr  3 14:12:29 2026 =-=-=-=-=
=-=-=-=-= Iteration 400 Fri Apr  3 14:12:31 2026 =-=-=-=-=
=-=-=-=-= Iteration 500 Fri Apr  3 14:12:33 2026 =-=-=-=-=
=-=-=-=-= Iteration 600 Fri Apr  3 14:12:35 2026 =-=-=-=-=
=-=-=-=-= Iteration 700 Fri Apr  3 14:12:37 2026 =-=-=-=-=
=-=-=-=-= Iteration 800 Fri Apr  3 14:12:40 2026 =-=-=-=-=
=-=-=-=-= Iteration 900 Fri Apr  3 14:12:42 2026 =-=-=-=-=

pred_r <- predict(model1, horizon = 12)
plot(pred_r)