bnlssm

Create Bayesian nonlinear non-Gaussian state-space model

Since R2023b

Description

bnlssm creates a bnlssm object, representing a Bayesian nonlinear non-Gaussian state-space model, from a specified nonlinear mapping function, which defines the state-space model structure, and the log prior distribution function of the parameters. The state-space model can be time-invariant or time-varying (see Decide on Model Structure), and the state or observation variables, x_t or y_t, respectively, can be multivariate series. State disturbances are Gaussian random variables, and the observation innovations have a custom distribution.

In general, a bnlssm object specifies the joint prior distribution and characteristics of the state-space model only. That is, the model object is a template intended for further use.

Alternative state-space models include:

• The ssm model object — Standard linear Gaussian state-space model

• The dssm model object — Standard linear Gaussian state-space model with diffuse initial state distribution

• The bssm model object — Bayesian linear state-space model

Creation

Syntax

PriorMdl = bnlssm(ParamMap,ParamDistribution)

PriorMdl = bnlssm(ParamMap,ParamDistribution,Name=Value)

Description

PriorMdl = bnlssm(ParamMap,ParamDistribution) creates the Bayesian nonlinear non-Gaussian state-space model object PriorMdl. ParamMap is a function of the collection of state-space model parameters Θ that characterizes the nonlinear state dynamics (transition of states x_t from time t – 1 to t) and nonlinear state measurements (observations) y_t. ParamDistribution is the log prior density of Θ. The state model has additive linear Gaussian disturbances u_t; the observation model can have additive linear Gaussian innovations ε_t or y_t can have a custom density.

PriorMdl is a template that specifies the joint prior distribution of Θ and the structure of the nonlinear state-space model.

PriorMdl = bnlssm(ParamMap,ParamDistribution,Name=Value) sets properties using name-value arguments. For example, bnlssm(ParamMap,ParamDistribution,ObservationForm="distribution") specifies that the observation log density log(p(y_t|x_t)) is custom and defined in ParamMap.

Input Arguments

ParamMap — Parameter Θ mapping characterizing state-space model structure

function handle

Parameter Θ mapping characterizing the state-space model structure and determining the data likelihood, specified as a function handle in the form @fcnName, where fcnName is the function name. ParamMap sets the ParamMap property.

This table contains the supported forms for ParamMap and their required signatures (For more details on the terms in the equations, see Bayesian nonlinear non-Gaussian state-space model).

The distribution of the observations y_t determines the form and paramMap is the MATLAB^® function name for ParamMap, but you can use a different name.

Observation Form	Observation Distribution	State-Space Model	Required ParamMap Signature
Equation	εt is an additive linear Gaussian random series.	$$\begin{array}{l}{x}_t=A_t\left(x_{t-1};{\theta }_A\right)+B_t\left({\theta }_B\right)u_t\\ y_t=C_t\left(x_t;{\theta }_C\right)+D_t\left({\theta }_D\right){\epsilon }_t.\end{array}$$	function [A,B,C,D,Mean0,Cov0] = paramMap(theta)
Distribution You must set ObservationForm="distribution"	Custom	$$\begin{array}{l}{x}_t=A_t\left(x_{t-1};{\theta }_{A_t}\right)+B_t\left({\theta }_{B_t}\right)u_t\\ y_t~\log p\left(y_t|x_t;{\theta }_{y_t}\right).\end{array}$$	function [A,B,LogY,Mean0,Cov0] = paramMap(theta)

paramMap can accept additional inputs, such as predictor data for a regression component in the observation equation, and it can return additional outputs. This signature shows the reserved, but optional, outputs.

function [A,B,LogY,Mean0,Cov0,StateType,DeflatedData] = paramMap(theta)

• theta is a numParams-by-1 numeric vector of the state-space model parameters Θ as the first input argument. The function accepts inputs in subsequent positions.

• ParamMap returns the state-space model parameters in this table.

  Equation Quantity	Output Position	Output
  State-transition mapping, At	1	Required output A is one of the following quantities: For a linear time-invariant model, A is an m-by-m coefficient matrix. For a linear time-varying model, A is a T-by-1 cell vector, where cell t contains the mt-by-mt–1 coefficient matrix. For a nonlinear mapping, A is a function handle. The corresponding function must have this signature: function xt = at(lagxt) lagxt is a column vector or matrix of states at time t – 1, with mt–1 rows. xt is a column vector or matrix at time t, with mt rows. The function can accept additional inputs and return additional outputs. For details on input and output matrices, see Multipoint.
  State-disturbance-loading coefficient matrix, Bt	2	Required output B is one of the following mappings for the additive linear Gaussian state-disturbance series ut: For a linear time-invariant model, B is an m-by-k coefficient matrix. For a linear time-varying model, B is a T-by-1 cell vector, where cell t contains the mt-by-kt coefficient matrix.
  Measurement-sensitivity mapping, Ct	3, for equation form only	Required output C, for equation-form models, is one of the following quantities: For a linear time-invariant model, C is an n-by-m coefficient matrix. For a linear time-varying model, C is a T-by-1 cell vector, where cell t contains the nt-by-m coefficient matrix. For a nonlinear mapping, C is a function handle. The corresponding function must have this signature: function yt = ct(xt) xt is a column vector or matrix of states at time t, with mt rows. yt is a column vector or matrix at time t, with mt rows. The function can accept additional inputs and return additional outputs. For details on input and output matrices, see Multipoint.
  Observation-innovation coefficient matrix, Dt	4, for equation form only	Required output D, for equation-form models, is one of the following mappings for the additive linear Gaussian observation-innovation series εt: For a linear time-invariant model, D is an n-by-h coefficient matrix. For a linear time-varying model, D is a T-by-1 cell vector, where cell t contains the nt-by-ht coefficient matrix.
  log p(yt|xt;Θyt)	3, for distribution form only	Required output LogY, for distribution-form models, is a function handle to the custom log observation density. The corresponding function must have this signature: function p = logyt(yt,xt) yt is a column vector of responses at time t yt with nt rows. xt is a column vector or matrix of states at time t, with mt rows. p is a scalar or row vector of corresponding log densities. The function can accept additional inputs and return additional outputs. For details on input and output matrices, see Multipoint.
  Initial state mean vector, μ0	5, for equation form 4, for distribution form	Output Mean0 is an m0-by-1 vector. For linear state transition A, the default is the mean of the stationary distribution of the states. For nonlinear state transitions, you must specify Mean0. bnlssm assumes x0 ~ N(μ0,Σ0) regardless of equation form. For more details, see State Characteristics.
  Initial state covariance matrix, Σ0	6	Output Cov0 is an m0-by-m0 matrix. For linear state transition A, the default is the covariance of the stationary distribution of the states. For nonlinear state transitions, you must specify Cov0. For more details, see State Characteristics.
  State classification vector, StateType	7	Optional output vector of flags specifying the corresponding state type, either stationary (0), the constant 1 (1), or diffuse, static, or nonstationary (2). For more details, see State Characteristics.
  DeflatedData	8	Optional output array of response data deflated by predictor data, which accommodates a regression component in the observation equation

  The subscript t of functions and parameters indicate that the parameters can be time-varying. Ignore the subscript for time-invariant functions or parameters.

  To skip specifying an optional output argument, set the argument to [] in the function body. For example, to skip specifying StateType, set StateType = []; in the function.

Specify parameters to include in the posterior distribution by setting their value to an entry in the first input argument theta and set known entries of the coefficients to their values. For example, the following lines define the 1-D time-invariant state-space model

$$\begin{array}{l}{x}_t=a{x}_{t-1}+b{u}_t\\ y_t=x_t+d{\epsilon }_t.\end{array}$$

A = theta(1);
B = theta(2);
C = 1;
D = theta(3);

If paramMap requires the input parameter vector argument only, you can create the bnlssm object by calling:

Mdl = bnlssm(@paramMap,...)

In general, create the bnlssm object by calling:

Mdl = bnlssm(@(theta)paramMap(theta,...otherInputArgs...),...)

Example: bnlssm(@(params)paramFun(theta,y,z),@ParamDistribution) specifies the function paramFun that accepts the state-space model parameters theta, observed responses y, and predictor data z.

Tip

• Because out-of-bounds prior density evaluation is 0, set the log prior density of out-of-bounds parameter arguments to -Inf.

• A best practice is to set StateType of each state within ParamMap for both of the following reasons:

  • By default, the software generates StateType, but the default choice might not be accurate. For example, the software cannot distinguish between a constant 1 state and a static state.

  • The software cannot infer StateType from data because the data theoretically comes from the observation equation. The realizations of the state equation are unobservable.

Data Types: function_handle

ParamDistribution — Log of joint probability density function of the state-space model parameters Π(Θ)

function handle

Log of joint probability density function of the state-space model parameters Π(Θ), specified as a function handle in the form @fcnName, where fcnName is the function name. ParamDistribution sets the ParamDistribution property.

Suppose logPrior is the name of the MATLAB function defining the joint prior distribution of Θ. Then, logPrior must have this form.

function logpdf = logPrior(theta,...otherInputs...)
	...
end

where:

• theta is a numparams-by-1 numeric vector of the linear state-space model parameters Θ. Elements of theta must correspond to those of ParamMap. The function can accept other inputs in subsequent positions.

• logpdf is a numeric scalar representing the log of the joint probability density of Θ at the input theta.

If ParamDistribution requires the input parameter vector argument only, you can create the bnlssm object by calling:

Mdl = bnlssm(...,@logPrior)

In general, create the bnlssm object by calling:

Mdl = bnlssm(...,@(theta)logPrior(theta,...otherInputArgs...))

Tip

Because out-of-bounds prior density evaluation is 0, set the log prior density of out-of-bounds parameter arguments to -Inf.

Data Types: function_handle

Properties

ObservationForm — Observation model form

"equation" (default) | "distribution"

Observation model form, specified as a value in this table.

Value	Observation Distribution	State-Space Model
"equation"	εt is an additive linear Gaussian random series.	$$\begin{array}{l}{x}_t=A_t\left(x_{t-1};{\theta }_A\right)+B_t\left({\theta }_B\right)u_t\\ y_t=C_t\left(x_t;{\theta }_C\right)+D_t\left({\theta }_D\right){\epsilon }_t.\end{array}$$
"distribution"	Custom	$$\begin{array}{l}{x}_t=A_t\left(x_{t-1};{\theta }_{A_t}\right)+B_t\left({\theta }_{B_t}\right)u_t\\ y_t~\log p\left(y_t|x_t;{\theta }_{y_t}\right).\end{array}$$

Example: ObservationForm="distribution"

Data Types: char | string

Multipoint — Multipoint evaluation of nonlinear functions

empty array ([]) (default) | "A" | "C" | "LogY" | string vector

Multipoint evaluation of nonlinear functions of ParamMap A, C and LogY, specified as a "A", "C", "LogY", or a string vector of such values. Specify Multipoint to speed up particle filtering routines.

For the specified nonlinear functions, bnlssm object functions can evaluate multiple points simultaneously. For example, suppose x₁ and x₂ are two points (state particles) to be evaluated by A_t(x_t) and Multipoint="A". The function At evaluates the concatenated points A_t([x₁ x₂]) = [Z₁ Z₂].

You must write the functions so that they can evaluate numparticles points, or states, simultaneously. To write nonlinear functions to support multipoint evaluation:

• For A = A_t, the function must accept an m_t–1-by-numparticles matrix of state particles at time t, and return a 1-by-numparticles row vector of corresponding evaluations. For example, A = x(1,:)./x(2,:).

• For C = C_t, the function must accept an n_t-by-1 column vector responses at time t and an m_t-by-numparticles matrix of state particles at time t, and return a 1-by-numparticles row vector of corresponding evaluations. At time t, the software applies each observation to all state particles. For example, C = theta(2)*x(1,:).*x(2,:).

• For LogY = log p(y_t|x_t;Θ_yt), the function must accept an n_t-by-1 column vector responses at time t and an m_t-by-numparticles matrix of state particles at time t, and return a 1-by-numparticles row vector of corresponding log density evaluations. At time t, the software applies each observation to all state particles.

If you disable multipoint evaluation (the default), functions process points sequentially, for example, functions evaluate A_t(x₁) = Z₁, and then they evaluate A_t(x₂) = Z₂.

When A_t or C_t are coefficient matrices, functions always apply multipoint evaluation because, for example, A_t[x₁ x₂] is well-defined.

Example: Multipoint=["A" "LogY"] indicates that bnlssm object functions can evaluate multiple points of the nonlinear functions A and LogY simultaneously.

Data Types: char | string

ParamMap — Parameter Θ mapping characterizing state-space model structure

function handle

Parameter Θ mapping characterizing state-space model structure, stored as a function handle and set by the ParamMap input argument. ParamMap completely specifies the structure of the state-space model.

Data Types: function_handle

ParamDistribution — Parameter distribution representation

function handle | numeric matrix

Parameter distribution representation, stored as a function handle or a numparams-by-numdraws numeric matrix.

• ParamDistribution is a function handle for the log prior distribution of the parameters ParamDistribution when you create PriorMdl directly by using bnlssm.

• ParamDistribution is a numparams-by-numdraws numeric matrix containing random draws from the posterior distribution of the parameters when you sample from the posterior using an object function. Rows correspond to the elements of theta and columns correspond to subsequent draws of the pseudo-marginal and particle-marginal Metropolis-Hastings samplers [1][2][3].

Data Types: function_handle

Object Functions

estimate	Estimate posterior distribution of Bayesian nonlinear non-Gaussian state-space model parameters
filter	Forward recursion of Bayesian nonlinear non-Gaussian state-space model
smooth	Backward recursion of Bayesian nonlinear non-Gaussian state-space model
simsmooth	Bayesian nonlinear non-Gaussian state-space model simulation smoother
simulate	Simulate posterior draws of Bayesian nonlinear non-Gaussian state-space model parameters
tune	Tune Bayesian nonlinear non-Gaussian state-space model posterior sampler

Examples

Create Bayesian Nonlinear Model in Equation Form

This example shows how to create the following Bayesian nonlinear state-space model in equation form by using bnlssm. The state-space model contains two independent, stationary, autoregressive states each with a model constant. The observations are a nonlinear function of the states with Gaussian noise. The prior distribution of the parameters is flat. Symbolically, the system of equations is

$$\left[\begin{array}{c}{x}_{t,1}\\ x_{t,2}\\ x_{t,3}\\ x_{t,4}\end{array}\right]=\left[\begin{array}{cccc}{\theta }_1& {\theta }_2& 0& 0\\ 0& 1& 0& 0\\ 0& 0& {\theta }_3& {\theta }_4\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{t-1,1}\\ x_{t-1,2}\\ x_{t-1,3}\\ x_{t-1,4}\end{array}\right]+\left[\begin{array}{cc}{\theta }_5& 0\\ 0& 0\\ 0& {\theta }_6\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_{t,1}\\ u_{t,3}\end{array}\right]$$

$$y_t=\log (\exp (x_{t,1}-{\mu }_1)+\exp (x_{t,3}-{\mu }_3))+{\theta }_7{\epsilon }_t.$$

${\mu }_1$ and ${\mu }_3$ are the unconditional means of the corresponding states. The initial distribution moments of each state are their unconditional mean and covariance.

Create a Bayesian nonlinear state-space model characterized by the system. The observation equation is in equation form, that is, the function composing the states is nonlinear and the innovation series ${\epsilon }_{\mathit{t}}$ is additive, linear, and Gaussian. The Local Functions section contains two functions required to specify the Bayesian nonlinear state-space model: the state-space model parameter mapping function and the prior distribution of the parameters. You can use the functions only within this script.

Mdl = bnlssm(@paramMap,@priorDistribution)

Mdl = 
  bnlssm with properties:

             ParamMap: @paramMap
    ParamDistribution: @priorDistribution
      ObservationForm: "equation"
           Multipoint: [1x0 string]

Mdl is a bnlssm model specifying the state-space model structure and prior distribution of the state-space model parameters. Because Mdl contains unknown values, it serves as a template for posterior analysis with observations.

Local Functions

These functions specify the state-space model parameter mappings, in equation form, and log prior distribution of the parameters.

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(theta)
    A = @(x)blkdiag([theta(1) theta(2); 0 1],[theta(3) theta(4); 0 1])*x; 
    B = [theta(5) 0; 0 0; 0 theta(6); 0 0];
    C = @(x)log(exp(x(1)-theta(2)/(1-theta(1))) + ...
        exp(x(3)-theta(4)/(1-theta(3))));
    D = theta(7);
    Mean0 = [theta(2)/(1-theta(1)); 1; theta(4)/(1-theta(3)); 1];         
    Cov0 = diag([theta(5)^2/(1-theta(1)^2) 0 theta(6)^2/(1-theta(3)^2) 0]);          
    StateType = [0; 1; 0; 1];     % Stationary state and constant 1 processes
end

function logprior = priorDistribution(theta)
    paramconstraints = [(abs(theta([1 3])) >= 1) (theta(5:7) <= 0)];
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

Create Bayesian Nonlinear Model in Distribution Form

This example shows how to create the following Bayesian nonlinear state-space model in distribution form by using bnlssm. The state-space model contains two independent, nonstationary states. The states combine to form the probability of success for a series of Bernoulli observations. The two states contain white Gaussian noise. The prior distribution of the parameters is flat. Symbolically, the system of equations is

$$\left[\begin{array}{c}{x}_{1,t}\\ x_{2,t}\end{array}\right]=\left[\begin{array}{cc}{\theta }_1& 0\\ 0& {\theta }_2\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\end{array}\right]+\left[\begin{array}{cc}{\theta }_3& 0\\ 0& {\theta }_4\end{array}\right]\left[\begin{array}{c}{u}_{1,t}\\ u_{2,t}\end{array}\right]$$

$$y_t\sim p(y_t;x_t).$$

$\mathit{p}\left({\mathit{y}}_{\mathit{t}};{\mathit{x}}_{\mathit{t}}\right)$ is the Bernoulli probability distribution with probability of success $\frac{\exp \left({\mathit{x}}_{1,\mathit{t}}\right)}{\exp \left({\mathit{x}}_{1,\mathit{t}}\right)+\exp \left({\mathit{x}}_{2,\mathit{t}}\right)}$. Arbitrarily assume that the initial distribution of each state has a mean and standard deviation of 0.5, and assume the states are initially uncorrelated.

Create a Bayesian nonlinear state-space model characterized by the system. The observation equation is in distribution form, that is, the function representing the observation model is a probability distribution. Specify that the observation model is in distribution form.

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

Mdl = bnlssm(@paramMap,@priorDistribution,ObservationForm="distribution")

Mdl = 
  bnlssm with properties:

             ParamMap: @paramMap
    ParamDistribution: @priorDistribution
      ObservationForm: "distribution"
           Multipoint: [1x0 string]

Local Functions

These functions specify the state-space model parameter mappings, in equation form, and log prior distribution of the parameters.

function [A,B,LogY,Mean0,Cov0,StateType] = paramMap(theta)
    A = @(x)[theta(1) 0; 0 theta(2)]*x; 
    B = [theta(3) 0; 0 theta(4)];
    LogY = @(y,x)sum(log(binopdf(y,1,exp(x(1))/(exp(x(2))+exp(x(1))))));
    Mean0 = [0.5; 0.5];
    Cov0 = eye(2);
    StateType = [2; 2];     % Nonstationary state processes
end

function logprior = priorDistribution(theta)
    paramconstraints = (theta(3:4) <= 0);
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

Prepare Model for Multipoint Evaluation

Multipoint evaluation can speed up particle filtering routines. This example shows how to prepare the Bayesian nonlinear state-space model in Create Bayesian Nonlinear Model in Distribution Form for multipoint evaluation.

Because A is a linear function, it is prepared for multipoint evaluation. However, the following line in the local functions of Create Bayesian Nonlinear Model in Distribution Form does not evaluate points columnwise with respect to input x.

...   
    LogY = @(y,x)sum(log(binopdf(y,1,exp(x(1))/(exp(x(2))+exp(x(1))))));
...

You can rewrite the log density of ${\mathit{y}}_{\mathit{t}}$ so that it evaluates the Bernoulli density for each state particle (column) of x.

...
    LogY = @(y,x)logbernpdf(y,x);
...

function bpdf = logbernpdf(y,x)
    p = exp(x(1,:))./(exp(x(2,:)) + exp(x(1,:)));    
    bpdf = sum((y.*log(p)+(1-y).*log(1-p)),1);
end

Create a Bayesian nonlinear state-space model characterized by the system.

• The observation equation is in distribution form, that is, the function representing the observation model is a probability distribution. Specify that the observation model is in distribution form.

• The function A and LogY are in a form that allows for simultaneous function evaluations during the particle filtering routines of bnlssm object functions.

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

Mdl = bnlssm(@paramMap,@priorDistribution,ObservationForm="distribution", ...
    Multipoint=["A" "LogY"])

Mdl = 
  bnlssm with properties:

             ParamMap: @paramMap
    ParamDistribution: @priorDistribution
      ObservationForm: "distribution"
           Multipoint: ["A"    "LogY"]

Local Functions

These functions specify the state-space model parameter mappings, in equation form, and log prior distribution of the parameters.

function [A,B,LogY,Mean0,Cov0,StateType] = paramMap(theta)
    A = @(x)[theta(1) 0; 0 theta(2)]*x; 
    B = [theta(3) 0; 0 theta(4)];
    LogY = @(y,x)logbernpdf(y,x);
    Mean0 = [1; 1];         
    Cov0 = eye(2);          
    StateType = [2; 2];     % Nonstationary state process
end

function bpdf = logbernpdf(y,x)
    p = exp(x(1,:))./(exp(x(2,:)) + exp(x(1,:)));
    bpdf = sum((y.*log(p)+(1-y).*log(1-p)),1);
end

function logprior = priorDistribution(theta)
    paramconstraints = (theta(3:4) <= 0);
    if(sum(paramconstraints))
        logprior = -Inf;
    else 
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

More About

Bayesian Nonlinear Non-Gaussian State-Space Model

A Bayesian nonlinear non-Gaussian state-space model is a Bayesian view of a state-space model with possibly nonlinear dynamics and non-Gaussian observations.

In general, a Bayesian nonlinear, multivariate, time-varying, non-Gaussian state-space model is the system of equations

$$\begin{array}{l}{x}_t=A_t\left(x_{t-1};\Theta \right)+B_t\left(\Theta \right)u_t\\ y_t~\log p\left(y_t|x_t;\Theta ,Z_t\right).\end{array}$$

for t = 1, ..., T and where:

• $$x_t=\left[x_{1,t},\mathrm{...},x_{m_t,t}\right]\prime$$ is an m_t-dimensional state vector describing the dynamics of some, possibly unobservable, phenomenon at period t. The initial state distribution x₀ ~ N(μ₀, Σ₀), where Mean0 = μ₀ and Cov0 = Σ₀.

• $$y_t=\left[y_{1,t},\mathrm{...},y_{n_t,t}\right]\prime$$ is an n_t-dimensional observation vector describing how the states are measured by observers at period t.

• $${\epsilon }_t=\left[{\epsilon }_{1,t},\mathrm{...},{\epsilon }_{h_t,t}\right]\prime$$ is an h_t-dimensional white-noise vector of observation innovations at period t. All innovations are multivariate Gaussian distributed.

• p(y_t|x_t;Θ,Z_t) characterizes the distribution of the observations at time t. p(y_t|x_t;Θ,Z_t) can take this general form or p(y_t|x_t;Θ,Z_t) = C_t(x_t;Θ) + D_tε_t.

• $$u_t=\left[u_{1,t},\mathrm{...},u_{k_t,t}\right]\prime$$ is a k_t-dimensional white-noise vector of state disturbances at period t. All disturbances are multivariate Gaussian distributed.

• ε_t and u_t are uncorrelated.

• A_t(x_t–1) is the state-transition mapping, a possibly nonlinear function of the previous state x_t–1 that characterizes the dynamics of the state-space model at time t.

• B_t is the linear state-disturbance-loading coefficient matrix at time t for the additive Gaussian state disturbances u_t.

• C_t(x_t) is the measurement-sensitivity mapping, a possibly nonlinear function of the current state x_t that characterizes how the latent states are observed at time t.

• D_t is the linear observation-innovation coefficient matrix at time t for the additive Gaussian observation innovations ε_t.

• The vector of model parameters Θ are treated as random variables with a joint prior distribution Π(Θ) = ParamDistribution

• For time-invariant state-space models,

  • $$Z_t=\left[\begin{array}{cccc}{z}_{t1}& z_{t2}& \cdots & z_{td}\end{array}\right]$$ is row t of a T-by-d matrix of predictors Z. Each column of Z corresponds to a predictor, and each successive row to a successive period. If the observations are multivariate, then all predictors deflate each observation.

  • β is a d-by-n matrix of regression coefficients for Z_t.

• The parameters in the functions or coefficients A_t, B_t, C_t, D_t, and β (when present) are model parameters arbitrarily collected in the vector Θ.

State-Transition Mapping A_t

The state-transition mapping A_t is a function that specifies how the states x_t are expected to transition from period t – 1 to t, for all t = 1,...,T. In other words, the expected state-transition equation at period t is E(x_t|x_t–1) = A_t(x_t–1). The corresponding output A of the ParamMap input can have the following forms:

• For linear time-invariant transitions, A is an m-by-m matrix, where m is the number of state variables.

• For linear time-varying transitions, A is a series of matrices represented by a T-by-1 cell vector, where A{t} contains an m_t-by-m_{t – 1} state-transition coefficient matrix. If the number of state variables changes from period t – 1 to t, m_t ≠ m_{t – 1}.an m-by-m matrix.

• For time-invariant or time-varying transitions, A is a function handle corresponding to the function the governs how states transition for each time t – 1 to t.

State-Disturbance-Loading Coefficient Matrix B_t

The state-disturbance-loading coefficient matrix B_t is a matrix or cell vector of matrices that specifies the additive error structure of the state disturbances u_t in the state-transition equation from period t – 1 to t, for all t = 1,...,T. In other words, the state-transition equation at period t is x_t = A_tx_t–1 + B_tu_t.

For time-invariant state-space models, the output B is an m-by-k matrix, where m is the number of state variables and k is the number of state disturbances. The quantity B*(B') is the state-disturbance covariance matrix for all periods.

For time-varying state-space models, B is a T-dimensional cell array, where B{t} contains an m_t-by-k_t state-disturbance-loading coefficient matrix. If the number of state variables or state disturbances changes at period t, the matrix dimensions between B{t-1} and B{t} vary. The quantity B{t}*(B{t}') is the state-disturbance covariance matrix for period t.

Measurement-Sensitivity Mapping C_t

The measurement-sensitivity mapping C_t is a function that specifies how the states x_t are expected to combine at period t to form the observations, y_t, for all t = 1,...,T. In other words, the expected observation equation at period t is E(y_t|x_t) = C_t(x_t). The corresponding output C of the ParamMap input can have the following forms:

• For linear time-invariant model, C is an n-by-m matrix, where n is the number of observation variables and m is the number of state variables.

• For linear time-varying model, C is a series of matrices represented by a T-by-1 cell vector, where C{t} contains an n_t-by-m_t measurement-sensitivity coefficient matrix. If the number of state or observation variables changes at period t, the matrix dimensions between C{t-1} and C{t} vary.

• For a time-invariant or time-varying model, C is a function handle corresponding to the function the governs how states combine to form observations at each time t.

Observation-Innovation Coefficient Matrix D_t

The observation-innovation coefficient matrix D_t is a matrix or cell vector of matrices that specifies the additive error structure of the observation innovations ε_t in the observation equation at period t, for all t = 1,...,T. In other words, the observation equation at period t is y_t = C_tx_t + D_tε_t.

For time-invariant state-space models, the output D is an n-by-h matrix, where n is the number of observation variables and h is the number of observation innovations. The quantity D*(D') is the observation-innovation covariance matrix for all periods.

For time-varying state-space models, the output D is a T-dimensional cell array, where D{t} contains an n_t-by-h_t matrix. If the number of observation variables or observation innovations changes at period t, then the matrix dimensions between D{t-1} and D{t} vary. The quantity D{t}*(D{t}') is the observation-innovation covariance matrix for period t.

State Characteristics

Other state characteristics include initial state moments and a description of the dynamic behavior of each state.

You can specify the state characteristics by including them as output arguments in ParamMap in their appropriate positions.

• Mean0 — Initial state mean μ₀, an m-by-1 numeric vector, where m is the number of states in x₁.

• Cov0 — Initial state covariance matrix Σ₀, an m-by-m positive semidefinite matrix.

• StateType — State dynamic behavior indicator, an m-by-1 numeric vector. This table summarizes the available types of initial state distributions.

  Value	State Dynamic Behavior Indicator
  0	Stationary (for example, ARMA models)
  1	The constant 1 (that is, the state is 1 with probability 1)
  2	Diffuse, nonstationary (for example, random walk model, seasonal linear time series), or static state

  For example, suppose that the state equation has two state variables: The first state variable is an AR(1) process, and the second state variable is a random walk. Specify the initial distribution types by setting StateType=[0; 2]; within the ParamMap function.

Static State

A static state does not change in value throughout the sample, that is, $$P\left(x_{t+1}=x_t\right)=1$$ for all t = 1,...,T.

Tips

• Load the data to the MATLAB workspace before creating the model.

• Create the parameter-to-matrix mapping function and log prior distribution function each as their own file.

• The equation-form model is a special case of the distribution-form model, with LogY = @(y,x) log(mvnpdf(y,C(x),D*D')).

• Avoid an arbitrary choice of the initial state distribution. bnlssm functions generate the initial particles from the specified initial state distribution, which impacts the performance of the nonlinear filter. If the initial state specification is bad enough, importance weights concentrate on a small number of particles in the first SMC iteration, which might produce unreasonable filtering results. This vulnerability of the nonlinear model behavior contrasts with the stability of the Kalman filter for the linear model, in which the initial state distribution usually has little impact on the filter because the prior is washed out as it processes data.

Algorithms

• When A_t is a coefficient matrix, For each state variable j, default values of Mean0 and Cov0 depend on StateType(j):

  • If StateType(j) = 0 (stationary state), bnlssm generates the initial value using the stationary distribution. If you provide all values in the coefficient matrices (that is, your model has no unknown parameters), bnlssm generates the initial values. Otherwise, the software generates the initial values during estimation.

  • If StateType(j) = 1 (constant state), Mean0(j) is 1 and Cov0(j) is 0.

  • If StateType(j) = 2 (nonstationary or diffuse state), Mean0(j) is 0 and Cov0(j) is 1e7.

• For static states that do not equal 1 throughout the sample, the software cannot assign a value to the degenerate, initial state distribution. Therefore, set the StateType of static states to 2. Consequently, the software treats static states as nonstationary and assigns the static state a diffuse initial distribution.

• bnlssm models do not store observed responses or predictor data. Supply the data wherever necessary using the appropriate input or name-value pair arguments.

• DeflateY is the deflated-observation data, which accommodates a regression component in the observation equation. For example, in this function, which has a linear regression component, Y is the vector of observed responses and Z is the vector of predictor data.

  function [A,B,C,D,Mean0,Cov0,StateType,DeflateY] = ParamFun(theta,Y,Z)
  	...
  	DeflateY = Y - params(9) - params(10)*Z;
  	...
  end

References

[1] Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. "Particle Markov Chain Monte Carlo Methods." Journal of the Royal Statistical Society Series B: Statistical Methodology 72 (June 2010): 269–342. https://doi.org/10.1111/j.1467-9868.2009.00736.x.

[2] Andrieu, Christophe, and Gareth O. Roberts. "The Pseudo-Marginal Approach for Efficient Monte Carlo Computations." Ann. Statist. 37 (April 2009): 697–725. https://dx.doi.org/10.1214/07-AOS574.

[3] Durbin, J, and Siem Jan Koopman. Time Series Analysis by State Space Methods. 2nd ed. Oxford: Oxford University Press, 2012.

Version History

Introduced in R2023b