Main Content

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, and the state or observation variables, xt or yt, 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

Description

example

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 xt from time t – 1 to t) and nonlinear state measurements (observations) yt. ParamDistribution is the log prior density of Θ. The state model has additive linear Gaussian disturbances ut; the observation model can have additive linear Gaussian innovations εt or yt 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.

example

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(yt|xt)) is custom and defined in ParamMap.

Input Arguments

expand all

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 yt determines the form and paramMap is the MATLAB® function name for ParamMap, but you can use a different name.

Observation FormObservation DistributionState-Space ModelRequired ParamMap Signature
Equationεt is an additive linear Gaussian random series.

xt=At(xt1;θAt)+Bt(θBt)utyt=Ct(xt;θCt)+Dt(θDt)εt.

function [A,B,C,D,Mean0,Cov0] = paramMap(theta)

Distribution

You must set ObservationForm="distribution"

Custom

xt=At(xt1;θAt)+Bt(θBt)utyt~logp(yt|xt;θyt).

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 QuantityOutput PositionOutput
    State-transition mapping, At1

    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, Bt2

    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, Ct3, 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, Dt4, 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|xtyt)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(μ00) regardless of equation form. For more details, see State Characteristics.

    Initial state covariance matrix, Σ06

    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, StateType7

    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.

    DeflatedData8

    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

xt=axt1+butyt=xt+dεt.

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

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

expand all

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

ValueObservation DistributionState-Space Model
"equation"εt is an additive linear Gaussian random series.

xt=At(xt1;θAt)+Bt(θBt)utyt=Ct(xt;θCt)+Dt(θDt)εt.

"distribution"Custom

xt=At(xt1;θAt)+Bt(θBt)utyt~logp(yt|xt;θyt).

Example: ObservationForm="distribution"

Data Types: char | string

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 x1 and x2 are two points (state particles) to be evaluated by At(xt) and Multipoint="A". The function At evaluates the concatenated points At([x1 x2]) = [Z1 Z2].

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

  • For A = At, the function must accept an mt–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 = Ct, the function must accept an nt-by-1 column vector responses at time t and an mt-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(yt|xtyt), the function must accept an nt-by-1 column vector responses at time t and an mt-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 At(x1) = Z1, and then they evaluate At(x2) = Z2.

When At or Ct are coefficient matrices, functions always apply multipoint evaluation because, for example, At[x1 x2] 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

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

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

filterForward recursion of Bayesian nonlinear non-Gaussian state-space model
smoothBackward recursion of Bayesian nonlinear non-Gaussian state-space model
simsmoothBayesian nonlinear non-Gaussian state-space model simulation smoother

Examples

collapse all

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

[xt,1xt,2xt,3xt,4]=[θ1θ200010000θ3θ40001][xt-1,1xt-1,2xt-1,3xt-1,4]+[θ50000θ600][ut,1ut,3]

yt=log(exp(xt,1-μ1)+exp(xt,3-μ3))+θ7εt.

μ1 and μ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 ε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

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

[xt,1xt,2]=[θ100θ2][xt-1,1xt-1,2]+[θ300θ4][ut,1ut,2]

ytp(yt;xt).

p(yt;xt) is the Bernoulli probability distribution with probability of success exp(xt,1)exp(xt,1)+exp(xt,2). 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

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 yt 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

expand all

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')).

Algorithms

  • When At 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

See Also

Objects