filter

R2026b

Forward recursion of Bayesian linear state-space model

Since R2026b

collapse all in page

Syntax

[X,logL] = filter(PriorMdl,Y,params)

[X,logL] = filter(PriorMdl,Y,params,Name=Value)

[X,logL,Output] = filter(___)

Description

[X,logL] = filter(PriorMdl,Y,params) returns the means X of the posterior state filtering distribution for each sampling time in the input response data Y, and the corresponding loglikelihood logL. To return the results, the function performs forward recursion, or filtering, of the prior Bayesian linear state-space model PriorMdl. filter computes the posterior distribution of the states conditioned on the state-space model parameters params. This computation does not incorporate parameter uncertainty.

filter obtains filtering distribution estimates of a model with Gaussian state disturbances and observation innovations, regardless of the values of the StateDistribution and ObservationDistribution properties of PriorMdl.

example

[X,logL] = filter(PriorMdl,Y,params,Name=Value) specifies additional options using one or more name-value arguments. For example, filter(PriorMdl,Y,params,Univariate=true,SquareRoot=true) specifies the univariate treatment of a multivariate state series, and use of the square root filter algorithm of the Kalman filter.

example

[X,logL,Output] = filter(___) additionally returns the output structure array Output using any of the input argument combinations in the previous syntaxes. Output contains the following information associated with each observation in the input data:

Loglikelihood value associated with the input data and parameters
State filter and forecast estimates with uncertainty
Observation forecast estimates with uncertainty
Kalman gain matrix
Flags indicating the data used by the software for filtering

example

Examples

collapse all

Compute Filter State Estimates and Loglikelihood

Open Live Script

Simulate observed responses from a known state-space model, estimate the posterior distribution of the parameters, and then obtain state estimates from the posterior filtering distribution.

Suppose the following state-space model is a data generating process (DGP).

$[\begin{array}{cccccccccccccccccccc} x_{t, 1} \\ x_{t, 2} \end{array}] = [\begin{array}{cccccccccccccccccccc} 0.5 & 0 \\ 0 & - 0.75 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{t - 1, 1} \\ x_{t - 1, 2} \end{array}] + [\begin{array}{cccccccccccccccccccc} 1 & 0 \\ 0 & 0.5 \end{array}] [\begin{array}{cccccccccccccccccccc} u_{t, 1} \\ u_{t, 2} \end{array}]$

$y_{t} = [\begin{array}{cccccccccccccccccccc} 1 & 1 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{t, 1} \\ x_{t, 2} \end{array}] .$

Create a standard state-space model object ssm that represents the DGP.

trueTheta = [0.5; -0.75; 1; 0.5];
A = [trueTheta(1) 0; 0 trueTheta(2)];
B = [trueTheta(3) 0; 0 trueTheta(4)];
C = [1 1];
DGP = ssm(A,B,C);

Simulate a response path from the DGP, and obtain filtered state estimates.

rng(1); % For reproducibility
y = simulate(DGP,200);
XFilter = filter(DGP,y);

Suppose the structure of the DGP is known, but the state parameters trueTheta are unknown, explicitly

$[\begin{array}{cccccccccccccccccccc} x_{t, 1} \\ x_{t, 2} \end{array}] = [\begin{array}{cccccccccccccccccccc} ϕ_{1} & 0 \\ 0 & ϕ_{2} \end{array}] [\begin{array}{cccccccccccccccccccc} x_{t - 1, 1} \\ x_{t - 1, 2} \end{array}] + [\begin{array}{cccccccccccccccccccc} σ_{1} & 0 \\ 0 & σ_{2} \end{array}] [\begin{array}{cccccccccccccccccccc} u_{t, 1} \\ u_{t, 2} \end{array}]$

$y_{t} = [\begin{array}{cccccccccccccccccccc} 1 & 1 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{t, 1} \\ x_{t, 2} \end{array}] .$

Consider a Bayesian state-space model that represents the model with unknown parameters. Assume that the prior distributions of $ϕ_{1}$ , $ϕ_{2}$ , $σ_{1}^{2}$ , and $σ_{2}^{2}$ are independent Gaussian with mean 0.5 and variance 1.

The Local Functions section contains two functions required to specify the Bayesian state-space model. You can use the functions only within this script.

The paramMap function accepts a vector of the unknown state-space model parameters and returns the following quantities:

A = $[\begin{array}{cc} ϕ_{1} & 0 \\ 0 & ϕ_{2} \end{array}]$ .
B = $[\begin{array}{cc} σ_{1} & 0 \\ 0 & σ_{2} \end{array}]$ .
C = $[\begin{array}{cc} 1 & 1 \end{array}]$ .
D = 0.
Mean0 and Cov0 are empty arrays [], which specify the defaults.
StateType = $[\begin{array}{cc} 0 & 0 \end{array}]$ , indicating that each state is stationary.

The priorDistribution function accepts the same vector of unknown parameters, but returns the log prior density of the parameters at their current values. The function specifies that parameter values outside the parameter space have a log prior density of -Inf.

Create the Bayesian state-space model by passing function handles for paramMap and priorDistribution directly to bssm.

PriorMdl = bssm(@paramMap,@priorDistribution)

PriorMdl = 
Mapping that defines a state-space model:
    @paramMap

Log density of parameter prior distribution:
    @priorDistribution

PriorMdl is a bssm object representing the Bayesian state-space model with unknown parameters.

Estimate the posterior distribution using estimate. Specify a random set of positive values in [0,1] to initialize the Markov chain Monte Carlo (MCMC) algorithm. Tune the sampler by applying a 500 burn-in period, thinning the sample by a factor of 30, and obtaining 2000 draws (after processing). Return parameter estimates from the posterior distribution.

numParams = 4;
theta0 = rand(numParams,1);
[PosteriorMdl,estParams] = estimate(PriorMdl,y,theta0,Univariate=true, ...
    BurnIn=500,NumDraws=2000,Thin=30);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
         Optimization and Tuning        
      | Params0  Optimized  ProposalStd 
----------------------------------------
 c(1) |  0.6968    0.4459      0.0798   
 c(2) |  0.7662   -0.8781      0.0483   
 c(3) |  0.3425    0.9633      0.0694   
 c(4) |  0.8459    0.3978      0.0726   
 
             Posterior Distributions            
      |   Mean     Std   Quantile05  Quantile95 
------------------------------------------------
 c(1) |  0.4495  0.0822     0.3135      0.5858  
 c(2) | -0.8561  0.0587    -0.9363     -0.7468  
 c(3) |  0.9645  0.0744     0.8448      1.0863  
 c(4) |  0.4333  0.0860     0.3086      0.5889  
Proposal acceptance rate = 38.85%

PosteriorMdl is a bssm object representing the posterior distribution.

Obtain state estimates from the posterior filtering distribution and return the data loglikelihood.

[XFilterHat,logL] = filter(PosteriorMdl,y,estParams);

Plot the state estimates with the true values.

figure
tiledlayout(2,1)
nexttile
plot(XFilter(:,1),"-",LineWidth=2)
hold on
plot(XFilterHat(:,1),"--",LineWidth=2)
hold off
title("x_1")
legend("True","Filtered")
nexttile
plot(XFilter(:,2),"-",LineWidth=2)
hold on
plot(XFilterHat(:,2),"--",LineWidth=2)
hold off
title("x_2")
legend("True","Filtered")
sgtitle("True State Values and Filtered State Estimates")

Figure contains 2 axes objects. Axes object 1 with title x indexOf 1 baseline contains 2 objects of type line. These objects represent True, Filtered. Axes object 2 with title x indexOf 2 baseline contains 2 objects of type line. These objects represent True, Filtered.

Local Functions

This example uses the following functions. paramMap is the parameter-to-matrix mapping function, and priorDistribution is the log prior distribution of the parameters.

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(theta)
    A = [theta(1) 0; 0 theta(2)];
    B = [theta(3) 0; 0 theta(4)];
    C = [1 1];
    D = 0;              % No observation noise
    Mean0 = [];         % MATLAB uses default initial state mean
    Cov0 = [];          % MATLAB uses default initial state covariances
    StateType = [0; 0]; % Two stationary states
end

function logprior = priorDistribution(theta)
    paramconstraints = [(abs(theta(1)) >= 1) (abs(theta(2)) >= 1) ...
        (theta(3) < 0) (theta(4) < 0)];
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
        mu0 = 0.5*ones(numel(theta),1); 
        sigma0 = 1;
        p = normpdf(theta,mu0,sigma0);
        logprior = sum(log(p));
    end
end

Obtain Filtered Coefficient Estimates of Regression Equation

Open Live Script

Model the coefficients of a linear regression as a random walk within a Bayesian state-space model. Track the evolution of the coefficients as the system processes observations.

The Bayesian state-space model is

$\begin{array}{l} y_{t} = β_{0, t} + β_{1, t} x_{1, t} + β_{2, t} x_{2, t} + u_{t} \\ β_{t} = β_{t - 1} + ε_{t}, \end{array}$

where:

$y_{t}$ is the US three-month T-bill rate.
$x_{1, t}$ is the US CPI-based inflation rate.
$x_{2, t}$ is the US M2 money supply growth rate.
$β_{t} = [\begin{array}{ccc} β_{0, t} & β_{1, t} & β_{2 t} \end{array}]'$ is the vector of coefficients.
$ε_{t}$ is a series of 3-D iid Gaussian disturbances. The variance of dimension $j$ is $σ_{j}^{2}$ and the prior distribution of $σ_{j}^{2}$ is $I G (2, 0.5)$ .
$u_{t}$ is an iid Gaussian innovation series. The variance is $σ_{4}^{2}$ and the prior distribution of $σ_{4}^{2}$ is $I G (2, 0.5)$ .

Load the US macroeconomic data Data_USEconModel.mat. The variable DataTimeTable is a timetable that contains the series in the regression model, among other series. Extract the response and predictor variables CPIAUCSL, M2SL, and TB3MS, and remove all observations that contain at least one missing value.

load Data_USEconModel
DTT = rmmissing(DataTimeTable(:,["CPIAUCSL" "M2SL" "TB3MS"]));

Plot the series separately.

figure
tiledlayout(2,2)
for j = 1:3
    nexttile
    plot(DTT.Time,DTT{:,j})
    title("Series: " + DTT.Properties.VariableNames{j})
    xlabel("Time")
end

Figure contains 3 axes objects. Axes object 1 with title Series: CPIAUCSL, xlabel Time contains an object of type line. Axes object 2 with title Series: M2SL, xlabel Time contains an object of type line. Axes object 3 with title Series: TB3MS, xlabel Time contains an object of type line.

Stabilize the series by applying the first difference to the three-month T-bill series, and by converting the predictor series to rates.

y = diff(DTT.TB3MS);
X = DTT{:,["CPIAUCSL" "M2SL"]};
X = price2ret(X);

The Local Functions section contains two functions required to specify the Bayesian state-space model. You can use the functions only within this script.

The paramMap function accepts a vector of the four variances in the state-space model, the predictor data X, and the sample size $T$ , and returns the following quantities:

$A = [\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}]$ .
$B = [\begin{array}{ccc} σ_{1} & 0 & 0 \\ 0 & σ_{2} & 0 \\ 0 & 0 & σ_{3} \end{array}]$ .
$C$ is a $T$ -by-1 cell vector, where cell $t$ is $x_{t}' β_{t}$ .
$D = σ_{4}$ .
Mean0 and Cov0 are empty arrays [], which specify the defaults.
StateType = $[\begin{array}{ccc} 2 & 2 & 2 \end{array}]'$ , indicating that each state is nonstationary.

Create the Bayesian state-space model by passing function handles for paramMap and priorDistribution directly to bssm. Compute the sample size $T$ , and set values for the shape and scale of the inverse gamma distribution to variables.

T = numel(y);
a = 2.0; 
b = 0.5; 
PriorMdl = bssm(@(theta)paramMap(theta,X,T),@(theta)priorDistribution(theta,a,b));

PriorMdl is a bssm object representing the Bayesian state-space model.

Fit the model to the data. Choose the initial state of the parameters randomly from $I G (2, 0.5)$ . Because the observation innovations are uncorrelated, set the univariate treatment of the multivariate model for computational efficiency. For numerical stability, apply the square root filter.

rng(1,"twister")
params0 = 1./gamrnd(a,b,4,1);
[PosteriorMdl,estParams] = estimate(PriorMdl,y,params0, ...
    Univariate=true,SquareRoot=true)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
         Optimization and Tuning        
      | Params0  Optimized  ProposalStd 
----------------------------------------
 c(1) |  0.5400    0.0705      0.0209   
 c(2) |  5.5331    0.1663      0.0957   
 c(3) |  3.3037    0.1663      0.0957   
 c(4) |  1.9393    1.0591      0.0537   
 
            Posterior Distributions            
      |  Mean     Std   Quantile05  Quantile95 
-----------------------------------------------
 c(1) | 0.0771  0.0198    0.0488      0.1114   
 c(2) | 0.2858  0.1441    0.1046      0.5602   
 c(3) | 0.3139  0.1816    0.1033      0.6640   
 c(4) | 1.0640  0.0577    0.9741      1.1538   
Proposal acceptance rate = 42.10%

PosteriorMdl = 
  bssm with properties:

                   ParamMap: @(theta)paramMap(theta,X,T)
          ParamDistribution: [4×1000 double]
          StateDistribution: [1×1 struct]
    ObservationDistribution: [1×1 struct]

            Posterior Distributions            
      |  Mean     Std   Quantile05  Quantile95 
-----------------------------------------------
 c(1) | 0.0771  0.0198    0.0488      0.1114   
 c(2) | 0.2858  0.1441    0.1046      0.5602   
 c(3) | 0.3139  0.1816    0.1033      0.6640   
 c(4) | 1.0640  0.0577    0.9741      1.1538   

State-space model type: ssm

State vector length: 3
Observation vector length: 1
State disturbance vector length: 3
Observation innovation vector length: 1
Sample size supported by model: 200

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equations:
x1(t) = x1(t-1) + (0.08)u1(t)
x2(t) = x2(t-1) + (0.29)u2(t)
x3(t) = x3(t-1) + (0.31)u3(t)

Observation equation of period 1:
y1(t) = x1(t) + (4.82e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 2:
y1(t) = x1(t) + (4.80e-03)x2(t) + (8.80e-03)x3(t) + (1.06)e1(t)

Observation equation of period 3:
y1(t) = x1(t) + (5.46e-03)x2(t) + (3.70e-03)x3(t) + (1.06)e1(t)

Observation equation of period 4:
y1(t) = x1(t) + (5.36e-03)x3(t) + (1.06)e1(t)

Observation equation of period 5:
y1(t) = x1(t) + (6.78e-03)x2(t) + (9.64e-03)x3(t) + (1.06)e1(t)

Observation equation of period 6:
y1(t) = x1(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 7:
y1(t) = x1(t) + (6.73e-03)x2(t) + (0.01)x3(t) + (1.06)e1(t)

Observation equation of period 8:
y1(t) = x1(t) + (1.01e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 9:
y1(t) = x1(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 10:
y1(t) = x1(t) + (4.68e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

... (Set 'Period' parameter to see intermediate equations)

Observation equation of period 191:
y1(t) = x1(t) + (2.46e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 192:
y1(t) = x1(t) + (0.01)x2(t) + (0.01)x3(t) + (1.06)e1(t)

Observation equation of period 193:
y1(t) = x1(t) + (8.12e-03)x2(t) + (0.01)x3(t) + (1.06)e1(t)

Observation equation of period 194:
y1(t) = x1(t) + (6.54e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 195:
y1(t) = x1(t) + (0.02)x2(t) + (0.01)x3(t) + (1.06)e1(t)

Observation equation of period 196:
y1(t) = x1(t) + (9.07e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 197:
y1(t) = x1(t) + (0.02)x2(t) + (6.44e-03)x3(t) + (1.06)e1(t)

Observation equation of period 198:
y1(t) = x1(t) + (7.53e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)

Observation equation of period 199:
y1(t) = x1(t) - (0.03)x2(t) + (0.04)x3(t) + (1.06)e1(t)

Observation equation of period 200:
y1(t) = x1(t) + (5.36e-03)x2(t) + (0.02)x3(t) + (1.06)e1(t)


Initial state distribution:

Initial state means
 x1  x2  x3 
  0   0   0 

Initial state covariance matrix
     x1     x2     x3    
 x1  1e+07  0      0     
 x2  0      1e+07  0     
 x3  0      0      1e+07 

State types
    x1       x2       x3   
 Diffuse  Diffuse  Diffuse

estParams = 4×1

    0.0771
    0.2858
    0.3139
    1.0640

Obtain posterior estimates from the filtering distribution of the regression coefficients. Apply the univariate treatment and the square root filter.

BetaFilter = filter(PosteriorMdl,y,estParams,Univariate=true,SquareRoot=true);

BetaFilter is a T-by-3 matrix containing the posterior estimates of the filtering distribution of the regression coefficients.

Display the evolution of the coefficient estimates. Plot the estimates separately.

figure
tiledlayout(2,2)
nexttile
plot(BetaFilter(:,1))
title("\beta_{0,t}")
nexttile
plot(BetaFilter(:,2))
title("\beta_{1,t}")
nexttile
plot(BetaFilter(:,3))
title("\beta_{2,t}")

Figure contains 3 axes objects. Axes object 1 with title beta indexOf 0 ,t baseline contains an object of type line. Axes object 2 with title beta indexOf 1 ,t baseline contains an object of type line. Axes object 3 with title beta indexOf 2 ,t baseline contains an object of type line.

Early filter state estimates appear unstable, but the estimates seem to settle as observations are filtered through the model. The following conditions explain the early instability:

The initial parameter values params0 are randomly drawn.
The initial state covariance is diffuse for nonstationary states by default.
Filter estimates for period $t$ are based on data from period 1 through $t$ . Therefore, filter estimates for early periods are informed by few observations.

You can remedy the early instability by setting the initial state mean Mean0 to the standard linear regression coefficient estimates, and the initial state covariance Cov0 to a matrix of zeros.

Compare the final estimates to the estimates from fitting a standard linear model.

Mdl = fitlm(X,y);
[Mdl.Coefficients.Estimate BetaFilter(end,:)']

ans = 3×2

   -0.0061   -0.0612
   20.2826   21.1997
  -12.4488  -14.4608

The estimates are in fair agreement.

Local Functions

This example uses the following functions. paramMap is the parameter-to-matrix mapping function, and priorDistribution is the log prior distribution of the parameters.

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(theta,Z,T)
    A = eye(3);
    B = diag(theta(1:3));
    C = cell(T,1);
    for t = 1:T
        C{t} = [1 Z(t,:)];
    end
    D = theta(4);
    Mean0 = [];            % MATLAB uses default initial state mean
    Cov0 = [];             % MATLAB uses default initial state covariances
    StateType = [2; 2; 2]; % Three nonstationary states
end

function logprior = priorDistribution(theta,a,b)
    paramconstraints = theta < 0;
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
    p = zeros(4,1);
    for j = 1:numel(p)
        p(j) = a*log(b) - gammaln(a) - (a+1)*log(theta(j)) - b./theta(j);
    end
        logprior = sum(p);
    end
end

Obtain Filtered State Estimate Variances

Open Live Script

Obtain state filtering posterior distribution covariances of the regression coefficients in the model from the example Obtain Filtered Coefficient Estimates of Regression Equation.

Load and preprocess the US macroeconomic data, and then create and estimate the Bayesian state-space model.

load Data_USEconModel
DTT = rmmissing(DataTimeTable(:,["CPIAUCSL" "M2SL" "TB3MS"]));
y = diff(DTT.TB3MS);
X = DTT{:,["CPIAUCSL" "M2SL"]};
X = price2ret(X)*100;
T = numel(y);
a = 2.0; 
b = 0.5; 
PriorMdl = bssm(@(theta)paramMap(theta,X,T),@(theta)priorDistribution(theta,a,b));

rng(100,"twister")
params0 = 1./gamrnd(a,b,4,1);
options = optimoptions("fminunc",Display="off");
[PosteriorMdl,estParams] = estimate(PriorMdl,y,params0, ...
    Univariate=true,SquareRoot=true,Display=false,Options=options);

Obtain posterior estimates from the filtering distribution of the regression coefficients. Apply the univariate treatment and the square root filter. Return the output structure that contains posterior estimates and variances, and extract the state filtering posterior distribution covariances into a cell vector.

[BetaFilter,~,out] = filter(PosteriorMdl,y,estParams,Univariate=true, ...
    SquareRoot=true);
BetaFilterCov = {out.FilteredStatesCov};

BetaFilterCov is a T-by-1 cell vector, where cell $t$ contains the 3-by-3 covariance matrix of the regression coefficients based on the state filtering posterior distribution.

Plot approximate 95% confidence ellipse projections for the filtered regression coefficient estimates of the final period in the sample by passing the filter estimates and t = T (200) to the local function plotConfEllipseProjBeta (see Local Functions). To see confidence ellipses for different values of t, use the control.

t = 200;
plotConfEllipseProjBeta(BetaFilter,BetaFilterCov,t)

$Figure contains 3 axes objects. Axes object 1 with title beta indexOf 0 baseline vs beta indexOf 1 baseline, xlabel \beta_0, ylabel \beta_1 contains 2 objects of type line. One or more of the lines displays its values using only markers Axes object 2 with title beta indexOf 0 baseline vs beta indexOf 2 baseline, xlabel \beta_0, ylabel \beta_2 contains 2 objects of type line. One or more of the lines displays its values using only markers Axes object 3 with title beta indexOf 1 baseline vs beta indexOf 2 baseline, xlabel \beta_1, ylabel \beta_2 contains 2 objects of type line. One or more of the lines displays its values using only markers$

Local Functions

This example uses the following functions. paramMap is the parameter-to-matrix mapping function, priorDistribution is the log prior distribution of the parameters, and plotConfEllipseProjBeta plots confidence ellipses of the coefficients.

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(theta,Z,T)
    A = eye(3);
    B = diag(theta(1:3));
    C = cell(T,1);
    for t = 1:T
        C{t} = [1 Z(t,:)];
    end
    D = theta(4);
    Mean0 = [];            % MATLAB uses default initial state mean
    Cov0 = [];             % MATLAB uses default initial state covariances
    StateType = [2; 2; 2]; % Three nonstationary states
end

function logprior = priorDistribution(theta,a,b)
    paramconstraints = theta < 0;
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
        p = zeros(4,1);
        for j = 1:numel(p)
            p(j) = a*log(b) - gammaln(a) - (a+1)*log(theta(j)) - b./theta(j);
        end
        logprior = sum(p);
    end
end

function plotConfEllipseProjBeta(est,covar,t)
    p = width(est);
    alpha = 0.05;
    chi2val = chi2inv(1-alpha,p); 
    theta = linspace(0,2*pi,100);
    unitCircle = [cos(theta); sin(theta)];
    mn = est(t,:).';
    Sigma = covar{t}; 
    pairs = nchoosek(1:p,2);
    numpairs = height(pairs);
    projs = cell(numpairs,1);
    for k = 1:numpairs
        idx = pairs(k,:);
        S2 = Sigma(idx,idx);
        S2 = (S2+S2')/2;
        [V,D] = eig(S2);
        D = diag(max(diag(D),0));
        A = V*sqrt(D)*sqrt(chi2val);
        pts = (A*unitCircle) + mn(idx);
        projs{k} = pts.';
    end
    betas = "\beta_" + string(0:(p-1));
    pairNames = betas(pairs);
    titles = join(pairNames, " vs ");
    tiledlayout("flow",Padding="compact",TileSpacing="compact");
    for k = 1:numpairs
        nexttile
        P = projs{k};
        plot(P(:,1),P(:,2),"-",LineWidth=1.5)
        hold on
        plot(mn(pairs(k,1)),mn(pairs(k,2)),"r.",MarkerSize=20)
        xlabel(betas(pairs(k,1)))
        ylabel(betas(pairs(k,2)))
        title(titles(k))
        axis equal
        grid on
        hold off
    end
    sgtitle("Ellipsoid Projections at Period " + t)
end

Input Arguments

collapse all

`PriorMdl` — Prior Bayesian linear state-space model
`bssm` object

Prior Bayesian linear state-space model, specified as a bssm object returned by bssm or ssm2bssm.

The function handles of the properties PriorMdl.ParamDistribution and PriorMdl.ParamMap determine the prior distribution and the data likelihood, respectively. filter evaluates PriorMdl.ParamMap at params before computing the state posterior distribution.

filter ignores the values of the properties StateDistribution and ObservationDistribution of PriorMdl, and assumes that the state disturbances and observation innovations are Gaussian.

`Y` — Observed response data
numeric matrix | cell vector of numeric vectors

Observed response data, from which filter forms the posterior distribution, specified as a numeric matrix or a cell vector of numeric vectors.

If PriorMdl is time invariant with respect to the observation equation, Y is a T-by-n matrix. Each row of the matrix corresponds to a period, and each column corresponds to a particular observation in the model. T is the sample size and n is the number of observations per period. The last row of Y contains the latest observations.
If PriorMdl is time varying with respect to the observation equation, Y is a T-by-1 cell vector. Y{t} contains an n_t-dimensional vector of observations for period t, where t = 1, ..., T. The corresponding dimensions of the coefficient matrices, outputs of PriorMdl.ParamMap, C{t}, and D{t} must be consistent with the matrix in Y{t} for all periods. The last cell of Y contains the latest observations.

NaN elements indicate missing observations. For details on how the Kalman filter accommodates missing observations, see Algorithms.

Data Types: double | cell

`params` — State-space model parameters Θ
numeric vector

State-space model parameters Θ used to evaluate the parameter mapping Mdl.ParamMap, specified as a numparams-by-1 numeric vector. Elements of params must correspond to the elements of the first input arguments of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

Usually, you pass either the posterior means of the parameters or a draw from their posterior distribution to params. The second output of estimate contains the posterior means of the parameters, and the first output of simulate contains draws from their posterior distribution.

Data Types: double

Name-Value Arguments

collapse all

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: filter(PriorMdl,Y,params,Univariate=true,SquareRoot=true) specifies the univariate treatment of a multivariate state series and use of the square root filter algorithm of the Kalman filter.

`Univariate` — Flag for univariate treatment of multivariate series
`false` (default) | `true`

Flag for the univariate treatment of a multivariate series, specified as false or true.

Value	Description
`true`	Applies the univariate treatment of a multivariate series, also known as sequential filtering
`false`	Does not apply sequential filtering

The univariate treatment can accelerate and improve numerical stability of the Kalman filter. However, all observation innovations must be uncorrelated. That is, D_tD_t' must be diagonal, where D_t (t = 1, ..., T) is the output coefficient matrix D of PriorMdl.ParamMap.

Example: Univariate=true

Data Types: logical

`SquareRoot` — Flag for square root filter method
`false` (default) | `true`

Flag for the square root filter method, specified as false or true.

Value	Description
`true`	Applies the square root filter method for the Kalman filter
`false`	Does not apply the square root filter method

If you think the eigenvalues of the filtered state or forecasted observation covariance matrices are close to zero, then specify SquareRoot=true. The square root filter is robust to numerical issues arising from the finite precision of calculations, but requires more computational resources.

Example: SquareRoot=true

Data Types: logical

`Tolerance` — Forecast uncertainty threshold
`0` (default) | nonnegative scalar

Forecast uncertainty threshold, specified as a nonnegative scalar.

If the forecast uncertainty for a particular observation is less than Tolerance during numerical estimation, then the software removes the uncertainty corresponding to the observation from the forecast covariance matrix before its inversion.

To overcome numerical obstacles during estimation, It is best practice to set Tolerance to a small number, such as 1e-15.

Example: Tolerance=1e-15

Data Types: double

Output Arguments

collapse all

`X` — Means of posterior state filtering distribution E(x_t|y₁,…,y_t,Θ)
numeric matrix | cell vector of numeric vectors

Means of the posterior state filtering distribution E(x_t|y₁,…,y_t,Θ), returned as a T-by-m numeric matrix or a T-by-1 cell vector of numeric vectors.

Each row corresponds to a time point in the sample. The last row contains the latest filtered states.

If PriorMdl is time invariant, filter returns a matrix. Each column corresponds to a state variable x_t.

If PriorMdl is time varying, X is a cell vector. Cell t contains a column vector of filtered state estimates with length m_t. Each column corresponds to a state variable.

`logL` — Loglikelihood function value
scalar

Loglikelihood function value, returned as a scalar.

Missing observations do not contribute to the loglikelihood.

`Output` — Filtering results by period
structure array

Filtering results by period, returned as a structure array.

Output is a T-by-1 structure, where element t corresponds to the filtering result at time t.

The following table describes the fields of Output when Univariate is false (the default), and identifies any changes to the fields when Univariate is true.

Field	Description	Estimate
`LogLikelihood`	Scalar loglikelihood objective function value	Not applicable
`FilteredStates`	m_t-by-1 vector of the filtered states	$E (x_{t} \| y_{1}, ..., y_{t}, Θ_{t})$
`FilteredStatesCov`	m_t-by-m_t variance-covariance matrix of the filtered states	$V a r (x_{t} \| y_{1}, ..., y_{t}, Θ_{t})$
`ForecastedStates`	m_t-by-1 vector of the state forecasts	$E (x_{t} \| y_{1}, ..., y_{t - 1}, Θ_{t})$
`ForecastedStatesCov`	m_t-by-m_t variance-covariance matrix of the state forecasts	$V a r (x_{t} \| y_{1}, ..., y_{t - 1}, Θ_{t})$
`ForecastedObs`	n_t-by-1 forecasted observation vector When `Univariate` is `true`, only the first elements between univariate and multivariate treatments of the problem are equal.	$E (y_{t} \| y_{1}, ..., y_{t - 1}, Θ_{t})$
`ForecastedObsCov`	n_t-by-n_t variance-covariance matrix of forecasted observations When `Univariate` is `true`, this field contains an n-by-1 vector of forecasted observation variances. The first element of this vector is equivalent to `ForecastedObsCov(1,1)` when `Univariate` is `false`. The rest of the elements are not necessarily equivalent to their corresponding values in `ForecastedObsCov` when `Univariate` is `false`.	$V a r (y_{t} \| y_{1}, ..., y_{t - 1}, Θ_{t})$
`KalmanGain`	m_t-by-h_t adjusted Kalman gain matrix When `Univariate` is `true`, might have different values for a multivariate treatment of the problem.	Not applicable
`DataUsed`	n_t-by-1 logical vector indicating whether the software filters using a particular observation. For example, if observation i at time t is a `NaN`, then element i in `DataUsed` at time t is `0`.	Not applicable

Algorithms

filter computes filtered state estimates under the following conditions:

filter evaluates the state-space model parameters at the input values in params. That is, the filtering state distribution is conditioned on the specified model parameter values.
The state-disturbance and observation-innovation distributions are Gaussian, regardless of the values of the properties PriorMdl.StateDistribution and PriorMdl.ObservationDistribution.

As a result, the computation of the means of the posterior state filtering distribution amounts to a forward recursion of the state-space model by using the standard Kalman filter, Therefore, this filtering method does not account for model parameter uncertainty (see the filter function of the ssm object). To work with more flexible models and incorporate model parameter uncertainty, use simsmooth.

filter accommodates missing data by not updating the filtered state estimates that correspond to missing observations. In other words, assume a missing observation at period t. The state forecast for period t based on the previous t – 1 observations is equivalent to the filtered state for period t.

Version History

Introduced in R2026b

filter

Syntax

Description

Examples

Compute Filter State Estimates and Loglikelihood

Obtain Filtered Coefficient Estimates of Regression Equation

Obtain Filtered State Estimate Variances

Input Arguments

`PriorMdl` — Prior Bayesian linear state-space model
`bssm` object

`Y` — Observed response data
numeric matrix | cell vector of numeric vectors

`params` — State-space model parameters Θ
numeric vector

Name-Value Arguments

`Univariate` — Flag for univariate treatment of multivariate series
`false` (default) | `true`

`SquareRoot` — Flag for square root filter method
`false` (default) | `true`

`Tolerance` — Forecast uncertainty threshold
`0` (default) | nonnegative scalar

Output Arguments

`X` — Means of posterior state filtering distribution E(x_t|y₁,…,y_t,Θ)
numeric matrix | cell vector of numeric vectors

`logL` — Loglikelihood function value
scalar

`Output` — Filtering results by period
structure array

Algorithms

Version History

See Also

Objects

Functions

filter

Syntax

Description

Examples

Compute Filter State Estimates and Loglikelihood

Obtain Filtered Coefficient Estimates of Regression Equation

Obtain Filtered State Estimate Variances

Input Arguments

PriorMdl — Prior Bayesian linear state-space model bssm object

Y — Observed response data numeric matrix | cell vector of numeric vectors

params — State-space model parameters Θ numeric vector

Name-Value Arguments

Univariate — Flag for univariate treatment of multivariate series false (default) | true

SquareRoot — Flag for square root filter method false (default) | true

Tolerance — Forecast uncertainty threshold 0 (default) | nonnegative scalar

Output Arguments

X — Means of posterior state filtering distribution E(xt|y1,…,yt,Θ) numeric matrix | cell vector of numeric vectors

logL — Loglikelihood function value scalar

Output — Filtering results by period structure array

Algorithms

Version History

See Also

Objects

Functions

`PriorMdl` — Prior Bayesian linear state-space model
`bssm` object

`Y` — Observed response data
numeric matrix | cell vector of numeric vectors

`params` — State-space model parameters Θ
numeric vector

`Univariate` — Flag for univariate treatment of multivariate series
`false` (default) | `true`

`SquareRoot` — Flag for square root filter method
`false` (default) | `true`

`Tolerance` — Forecast uncertainty threshold
`0` (default) | nonnegative scalar

`X` — Means of posterior state filtering distribution E(x_t|y₁,…,y_t,Θ)
numeric matrix | cell vector of numeric vectors

`logL` — Loglikelihood function value
scalar

`Output` — Filtering results by period
structure array