Description mathématique du modèle SV-SPM

Le modèle de volatilité stochastique avec prior spatial (SV-SPM) intègre une structure spatiale dans la modélisation de la volatilité des séries financières. Il est particulièrement utile pour capturer les dépendances temporellles et spatiales dans les rendements d'actifs.

Équation d’observation

Fichier principle sv_spm_bayes_main.m

% Simulation et estimation MCMC d'un modèle SV-SPM
clear; clc; close all;

%% 1. Génération des données simulées
fprintf('Génération des données...\n');
[N, returns, params_true] = simulate_sv_spatial_data(800);

fprintf('Taille de l’échantillon : %d\n', N);
disp('Valeurs réelles des paramètres :');
fprintf('  α = %.4f\n', params_true.alpha);
fprintf('  ν = %.4f\n', params_true.nu);
fprintf('  ρ = %.4f\n', params_true.rho);
fprintf('  τ = %.4f\n', sqrt(params_true.tau_sq));
fprintf('  σ_g = %.4f\n', sqrt(params_true.sigma_g_sq));
fprintf('  ℓ = %.4f\n', params_true.length_scale);

%% 2. Configuration de la chaîne MCMC
config.iter_max = 6000;
config.burn_in = 1500;
config.thin_by = 5;
config.samples_saved = floor((config.iter_max - config.burn_in)/config.thin_by);

% Priors informatifs
priors.alpha_mean = 0;     priors.alpha_var = 100;
priors.nu_mean   = 0;      priors.nu_var   = 100;
priors.rho_a     = 15;     priors.rho_b    = 1.2;
priors.tau_shape = 2.0;    priors.tau_rate = 0.02;
priors.gp_shape  = 2.0;    priors.gp_rate  = 0.02;
priors.len_shape = 1.8;    priors.len_rate = 0.08;

% Variances des propositions (marche aléatoire)
props.alpha_var = 0.015;
props.nu_var    = 0.015;
props.rho_var   = 0.012;
props.tau_var   = 0.018;
props.gp_var    = 0.020;
props.len_var   = 0.010;

fprintf('Configuration MCMC :\n');
fprintf('  Itérations totales : %d\n', config.iter_max);
fprintf('  Phase de burn-in : %d\n', config.burn_in);
fprintf('  Échantillons conservés : %d\n', config.samples_saved);

%% 3. Initialisation des paramètres
psi.alpha     = median(returns);
psi.nu        = log(var(returns));
psi.rho       = 0.92;
psi.tau_sq    = 0.12;
psi.sigma_g_sq = 0.8;
psi.length_scale = 1.1;

s_state = psi.nu * ones(N,1);  % État initial de la volatilité

% Préallocation des chaînes
chains = struct();
chains.alpha     = zeros(config.samples_saved, 1);
chains.nu        = zeros(config.samples_saved, 1);
chains.rho       = zeros(config.samples_saved, 1);
chains.tau_sq    = zeros(config.samples_saved, 1);
chains.sigma_g_sq = zeros(config.samples_saved, 1);
chains.length_scale = zeros(config.samples_saved, 1);
chains.posterior = zeros(config.samples_saved, 1);

accept_rates = zeros(6,1);

%% 4. Exécution de l'échantillonnage MCMC
fprintf('Démarrage de l’échantillonnage...\n');
current_psi = psi;
latent_vol = s_state;
idx_save = 0;

for iter = 1:config.iter_max
    if mod(iter, 500) == 0
        fprintf('  Itération %d/%d, taux d’acceptation moyen : %.1f%%\n', ...
                iter, config.iter_max, mean(accept_rates)*100);
    end

    % Mise à jour de chaque composante
    [current_psi.alpha, acc] = update_alpha(returns, latent_vol, current_psi, priors, props);
    accept_rates(1) = accept_rates(1) + acc;

    [current_psi.nu, latent_vol, acc] = update_nu_with_gp(returns, latent_vol, current_psi, priors, props, N);
    accept_rates(2) = accept_rates(2) + acc;

    [current_psi.rho, acc] = update_rho(latent_vol, current_psi, priors, props);
    accept_rates(3) = accept_rates(3) + acc;

    [current_psi.tau_sq, acc] = update_tau_sq(latent_vol, current_psi, priors, props);
    accept_rates(4) = accept_rates(4) + acc;

    [current_psi.sigma_g_sq, acc] = update_gp_variance(current_psi.nu, priors, props);
    accept_rates(5) = accept_rates(5) + acc;

    [current_psi.length_scale, acc] = update_length_scale(current_psi.nu, priors, props);
    accept_rates(6) = accept_rates(6) + acc;

    % Rééchantillonnage de l'état caché via FFBS
    latent_vol = sample_latent_volatility(returns, current_psi, latent_vol, N);

    % Stockage post-burn-in
    if iter > config.burn_in && mod(iter, config.thin_by) == 0
        idx_save = idx_save + 1;
        chains.alpha(idx_save)        = current_psi.alpha;
        chains.nu(idx_save)           = current_psi.nu;
        chains.rho(idx_save)          = current_psi.rho;
        chains.tau_sq(idx_save)       = current_psi.tau_sq;
        chains.sigma_g_sq(idx_save)   = current_psi.sigma_g_sq;
        chains.length_scale(idx_save) = current_psi.length_scale;
        chains.posterior(idx_save)    = compute_log_posterior(returns, latent_vol, current_psi, priors, N);
    end
end

fprintf('Estimation terminée.\n');

%% 5. Analyse des résultats
analyze_mcmc_output(chains, params_true, config, returns, latent_vol);

%% 6. Sauvegarde
save('sv_spm_estimation_results.mat', 'chains', 'params_true', 'config', 'returns', 'latent_vol');
fprintf('Résultats sauvegardés.\n');

Générateur de données simulate_sv_spatial_data.m

function [T, r, theta] = simulate_sv_spatial_data(T)
theta.alpha = 0.03;
theta.nu = -0.6;
theta.rho = 0.97;
theta.tau_sq = 0.06;
theta.sigma_g_sq = 0.6;
theta.length_scale = 0.9;

positions = (1:T)';
D = squareform(pdist(positions));
Sigma_nu = theta.sigma_g_sq * exp(-D.^2 ./ (2*theta.length_scale^2)) + 1e-6*eye(T);
nu_vec = chol(Sigma_nu,'lower') * randn(T,1);

s = zeros(T,1);
s(1) = nu_vec(1);
for t = 2:T
    s(t) = nu_vec(t) + theta.rho*(s(t-1)-nu_vec(t-1)) + sqrt(theta.tau_sq)*randn();
end

r = theta.alpha + exp(s/2).*randn(T,1);
fprintf('Données générées : T=%d\n', T);
end

Fonctions d’échantillonnage Metropolis-Hastings

% Mise à jour de alpha
function [a_new, accepted] = update_alpha(y, s, psi, pr, prop)
a_curr = psi.alpha;
a_prop = a_curr + sqrt(prop.alpha_var)*randn();
log_ratio = log_likelihood_alpha(y, s, a_prop) - log_likelihood_alpha(y, s, a_curr) + ...
            log_prior_alpha(a_prop, pr) - log_prior_alpha(a_curr, pr);
accepted = log(rand()) < log_ratio;
a_new = accepted ? a_prop : a_curr;
end

% Mise à jour de nu avec GP
function [nu_new, s_new, accepted] = update_nu_with_gp(y, s, psi, pr, prop, T)
nu_curr = psi.nu;
nu_prop = nu_curr + sqrt(prop.nu_var)*randn();

psi_temp = psi; psi_temp.nu = nu_prop;
log_ratio = log_likelihood_nu_gp(y, s, psi_temp, T) - log_likelihood_nu_gp(y, s, psi, T) + ...
            log_prior_nu(nu_prop, pr) - log_prior_nu(nu_curr, pr);
accepted = log(rand()) < log_ratio;

if accepted
    nu_new = nu_prop;
    s_new = sample_latent_volatility(y, psi_temp, s, T);
else
    nu_new = nu_curr;
    s_new = s;
end
end

% Mise à jour de rho
function [r_new, accepted] = update_rho(s, psi, pr, prop)
r_curr = psi.rho;
r_prop = max(min(r_curr + sqrt(prop.rho_var)*randn(), 0.995), -0.995);
log_ratio = log_likelihood_rho(s, psi, r_prop) - log_likelihood_rho(s, psi, r_curr) + ...
            log_prior_rho(r_prop, pr) - log_prior_rho(r_curr, pr);
accepted = log(rand()) < log_ratio;
r_new = accepted ? r_prop : r_curr;
end

% Mise à jour de tau_sq
function [t_new, accepted] = update_tau_sq(s, psi, pr, prop)
t_curr = psi.tau_sq;
log_t_prop = log(t_curr) + sqrt(prop.tau_var)*randn();
t_prop = exp(log_t_prop);
log_ratio = log_likelihood_tau(s, psi, t_prop) - log_likelihood_tau(s, psi, t_curr) + ...
            log_prior_tau(t_prop, pr) - log_prior_tau(t_curr, pr);
accepted = log(rand()) < log_ratio;
t_new = accepted ? t_prop : t_curr;
end

Rééchantillonnage de la volatilité latente (FFBS)

function s_out = sample_latent_volatility(y, psi, s_in, T)
s_out = zeros(T,1);
mu_filt = zeros(T,1); var_filt = zeros(T,1);

% Filtrage avant
mu_filt(1) = psi.nu;
var_filt(1) = psi.tau_sq;

for t = 2:T
    mu_pred = psi.nu + psi.rho*(s_in(t-1) - psi.nu);
    var_pred = psi.rho^2 * var_filt(t-1) + psi.tau_sq;
    
    innov = y(t) - psi.alpha;
    F = exp(s_in(t));
    gain = var_pred / (var_pred + F);
    mu_filt(t) = mu_pred + gain*(innov - mu_pred);
    var_filt(t) = (1 - gain)*var_pred;
end

% Échantillonnage arrière
s_out(T) = mu_filt(T) + sqrt(var_filt(T))*randn();
for t = T-1:-1:1
    J = psi.rho * var_filt(t) / (psi.rho^2 * var_filt(t) + psi.tau_sq);
    s_out(t) = mu_filt(t) + J*(s_out(t+1) - psi.nu - psi.rho*(s_out(t+1) - psi.nu));
end
end

Fonctions de vraisemblance et de prior

function ll = log_likelihood_alpha(y, s, a)
ll = -0.5 * sum((y - a).^2 ./ exp(s));
end

function lp = log_prior_alpha(a, pr)
lp = -0.5*(a - pr.alpha_mean)^2 / pr.alpha_var;
end

function lp = log_prior_rho(r, pr)
lp = (pr.rho_a - 1)*log(r + 1) + (pr.rho_b - 1)*log(1 - r);
end

function lp = log_prior_tau(t, pr)
lp = (pr.tau_shape - 1)*log(t) - pr.tau_rate / t;
end

Analyse des sorties MCMC

function analyze_mcmc_output(chains, truth, cfg, data, vol)
figure('Color','w','Position',[100 100 1300 900]);

% Tracés des chaînes
subplot(3,3,1); plot(chains.alpha, 'k'); hold on; yline(truth.alpha, 'r--');
title('Chaîne pour α'); xlabel('Itération'); ylabel('Valeur'); grid on;

subplot(3,3,2); plot(chains.rho, 'k'); hold on; yline(truth.rho, 'r--');
title('Chaîne pour ρ'); xlabel('Itération'); ylabel('Valeur'); grid on;

subplot(3,3,4); histogram(chains.alpha, 40, 'Normalization','pdf'); hold on;
xline(truth.alpha, 'r--'); title('Distribution a posteriori de α');

subplot(3,3,5); histogram(chains.rho, 40, 'Normalization','pdf'); hold on;
xline(truth.rho, 'r--'); title('Distribution a posteriori de ρ');

subplot(3,3,7); plot(exp(vol/2), 'b'); title('Volatilité estimée');
xlabel('Temps'); ylabel('Volatilité'); grid on;

subplot(3,3,8); autocorr(chains.alpha, 40); title('ACF de α');

sgtitle('Analyse bayésienne du modèle SV-SPM', 'FontSize',14);

% Résumé statistique
fprintf('\n=== Résultats de l’estimation ===\n');
fields = {'alpha','nu','rho','tau_sq','sigma_g_sq','length_scale'};
names = {'α','ν','ρ','τ²','σ_g²','ℓ'};

for k = 1:length(fields)
    chain_k = chains.(fields{k});
    m = mean(chain_k); sd = std(chain_k);
    ci = prctile(chain_k, [2.5 97.5]);
    true_val = truth.(fields{k});
    fprintf('%s : %.4f ± %.4f [%.4f, %.4f], vrai=%.4f\n', ...
            names{k}, m, sd, ci(1), ci(2), true_val);
end
end