How do I plot 14th and 86th percentile?

Here I am currently running Bayesian time varying model in matlab and this is my whole code. I am not sure how can I create 14th and 86th percentile of my posterior distribution/credible interval that I have found using gibbs sampling and carter kohn algorithm. Furthermore, how do I plot the 14th and 86th percentile along with my mean posterior distribution in one graph?

clear;

% Load data from Excel file
Cyclicalitydata = readtable("Could be final data set.xlsx");

Error using readtable (line 517)
Unable to find or open 'Could be final data set.xlsx'. Check the path and filename or file permissions.

Countrydata = Cyclicalitydata(strcmp(Cyclicalitydata.Country, 'Netherlands'), :);

% Prepare variables
Fiscalvariable = Countrydata.CAPB1; % Dependent variable, original series
LaggedFiscalvariable = [NaN; Fiscalvariable(1:end-1)]; % Creating one-period lag
Economicactivity = Countrydata.Outputgap; % Independent main variable
Govsize = Countrydata.Govsize; % Control variable 1
Debtp = Countrydata.DebtGDP; % Control variable 2
CPIlevel = Countrydata.CPIlevel; % Additional control variable 3
Years = Countrydata.Year; % Year for plotting and analysis
CPIgrowth = [NaN; CPIlevel(1:end-1)];

% Define the state-space model components
T = numel(LaggedFiscalvariable); % Number of observations, using the lagged variable now
C = [Economicactivity, Govsize, Debtp, CPIgrowth]; % Each row of C corresponds to an observation time
N = size(C,2);
A = eye(N); % Assuming independent evolution of each coefficient
Q = 0.01 * eye(N); % Small variance to assume minor changes over time
R = 1; % Observation noise variance

% Calculate OLS estimates for initialization
validIdx = ~isnan(LaggedFiscalvariable) & ~any(isnan(C), 2);
C_valid = C(validIdx, :);
Y_valid = LaggedFiscalvariable(validIdx);
beta_OLS = (C_valid' * C_valid) \ (C_valid' * Y_valid);
residuals = Y_valid - C_valid * beta_OLS;
sigma2_OLS = residuals' * residuals / (length(Y_valid) - size(C_valid, 2)); 

% Initialize beta_0 and V_0 with OLS estimates
beta_0 = beta_OLS;  % Use OLS estimates as initial states
V_0 = sigma2_OLS * eye(N);  % Scale identity matrix by OLS residual variance

% Number of Gibbs sampling iterations
numIterations = 1000;
beta_samples = zeros(N, T, numIterations); % Store samples of all beta coefficients

% Gibbs Sampling Loop
for iter = 1:numIterations
    % Run the Carter-Kohn algorithm for the current iteration
    [beta_t, ~] = carterKohn(LaggedFiscalvariable, A, C, Q, R, beta_0, V_0, T);

    % Store the sample
    beta_samples(:, :, iter) = beta_t;
end

% Post-processing after Gibbs sampling
% Focus on the OutputGap coefficient's estimates and its posterior distribution
burn_in = 100;
posterior_beta_outputgap = squeeze(mean(beta_samples(1, :, burn_in:end), 3));

% Plot only the OutputGap coefficient estimates over years
figure;
plot(Years, posterior_beta_outputgap, 'LineWidth', 2);
title('Posterior Estimates of OutputGap Coefficient over Years for the Netherlands');
xlabel('Year');
ylabel('Coefficient for OutputGap');


function [beta_t, V_t] = carterKohn(Y, A, C, Q, R, beta_0, V_0, T)
    n = size(A, 1); % Number of coefficients
    beta_t = zeros(n, T);
    V_t = zeros(n, n, T);
    beta_pred = beta_0;
    V_pred = V_0;
    
    for t = 1:T
        if isnan(Y(t))
            continue; % Skip iterations where Y is NaN
        end
        
        % Observation update
        C_t = C(t, :);  % Current row of C
        y_pred = C_t * beta_pred;
        S = C_t * V_pred * C_t' + R;
        K = V_pred * C_t' / S;
        beta_t(:, t) = beta_pred + K * (Y(t) - y_pred);
        V_t(:, :, t) = (eye(n) - K * C_t) * V_pred;
        
        % Time update
        if t < T
            beta_pred = A * beta_t(:, t);
            V_pred = A * V_t(:, :, t) * A' + Q;
        end
    end
end

Answer 1

[Star Strider]

Use the prctile function to calculate those percentiles.

[Yohan]

Hi, could you maybe help me how can I use it exactly? Because I have tried using it but for some reason i get the same error where my credible interval for those percentiles are equal to my mean. Since I am new to this program and I used chatgpt for many of the code so I think I must be doing something wrong here.

[Star Strider]

I cannot demonstrate it on your data because they are missing.

An example would be —

v = rand(1,20)

v = 1x20

0.5781 0.4799 0.1543 0.1357 0.6852 0.4078 0.5822 0.4300 0.9955 0.2118 0.6767 0.4732 0.6517 0.0304 0.7683 0.8073 0.7897 0.0232 0.3163 0.0813

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

vp = prctile(v, [14 86])

vp = 1x2

0.0976 0.7833

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot([1;1].*v, [1;1].*ones(size(v)), 'sb')
xline(vp,'-r')
grid

With a 20-element vector, there should be 2.8 elements below the 14th percentile line, and 2.8 elements above the 86th percentile line. The plotted percentiles are acceptably close, and will be increasingly more accurate as the number of elements increases. Note that ‘v’ does not have to be sorted first.

.