Main Content

bssm

Create Bayesian state-space model

Since R2022a

Description

bssm creates a bssm object, representing a Bayesian linear state-space model, from a specified parameter-to-matrix 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.

In general, a bssm 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. Specifically, to incorporate data into the model for posterior distribution analysis, pass the model object and data to the appropriate object function.

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 bnlssm model object — Bayesian nonlinear non-Gaussian state-space model

Creation

Description

There are several ways to create a bssm object representing a Bayesian state-space model:

  • Create Prior Model — Directly create a bssm object, representing a prior model, by using the bssm function and specifying the parameter-to-matrix mapping function and log prior distribution function of the parameters. This method accommodates simple through complex state-space model structures. For more details, see the following syntax.

  • Convert Standard Model to Bayesian Prior Model — Convert a specified standard, linear state-space model, an ssm object, to a bssm object by passing the model to the ssm2bssm function. The converted Bayesian model has the same state-space structure as the standard model. Use this method for simple state-space models when you prefer to specify the coefficient matrices explicitly. You can optionally specify the log prior density of the parameters. For details, see ssm2bssm and ssm.

  • Estimate Posterior Model — Pass a bssm object representing a prior model, observed response data, and initial parameter values to the estimate function to obtain a bssm object representing a posterior model. For more details, see estimate.

example

PriorMdl = bssm(ParamMap,ParamDistribution) creates the Bayesian state-space model object PriorMdl, representing a prior model with Gaussian state disturbances and observation innovations, using the parameter-to-matrix mapping function ParamMap, which you write, and the log prior distribution of the parameters.

The ParamMap function maps the collection of linear state-space model parameters Θ to the time-invariant or time-varying coefficient matrices A, B, C, and D. In addition to the coefficient matrices, ParamMap can optionally map any of the following quantities:

  • Initial state means and covariances Mean0 and Cov0

  • State types StateType

  • Deflated observation data DeflatedData to accommodate a regression component in the observation equation

The ParamDistribution function accepts Θ and returns the corresponding log density.

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

example

PriorMdl = bssm(ParamMap,ParamDistribution,Name=Value) sets properties representing the distributions of the state disturbances and observation innovations using name-value arguments. For example, bssm(ParamMap,ParamDistribution,ObservationDistribution=struct("Name","t","DoF",6)) specifies t distribution with 6 degrees of freedom for all observation innovation variables ut in the Bayesian state-space model.

Input Arguments

expand all

Parameter-to-matrix mapping function that determines the data likelihood, specified as a function handle in the form @fcnName, where fcnName is the function name. ParamMap sets the ParamMap property. Object functions of bssm compute the data likelihood by using the standard Kalman filter, where probabilities additionally condition on the parameters.

Suppose paramMap is the name of the MATLAB® function specifying how the state-space model parameters Θ map to the coefficient matrices and, optionally, other state-space model characteristics. Then, paramMap must have this form.

function [A,B,C,D,Mean0,Cov0,StateType,DeflateY] = paramMap(theta,...otherInputs...)
	...
end
where:

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 bssm object by calling:

Mdl = bssm(@paramMap,...)

In general, create the bssm object by calling:

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

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

Tip

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 bssm object by calling:

Mdl = bssm(...,@logPrior)

In general, create the bssm object by calling:

Mdl = bssm(...,@(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

Parameter-to-matrix mapping function, 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 bssm or when you convert a standard state-space model by using ssm2bssm.

  • ParamDistribution is a numParams-by-numDraws numeric matrix containing random draws from the posterior distribution of the parameters when you estimate the posterior by using estimate. Rows correspond to the elements of theta and columns correspond to subsequent draws of the Markov chain Monte Carlo sampler, such as the Metropolis-Hastings sampler [1][2].

Data Types: function_handle

Since R2022b

Distribution of the state disturbance process, specified as a distribution name or structure array in this table. bssm stores the value as a structure array.

DistributionNameVariable SupportStructure ArrayHyperparameter Estimation Support
Standard Gaussian"Gaussian"Multivariatestruct("Name","Gaussian")Not applicable
Student’s t"t"Multivariatestruct("Name","t","DoF",dof)Yes

  • The specified distribution applies to all state disturbance process variables.

  • For the Student's t distribution:

    • If you supply a structure array, you must specify both the "Name" and "DoF" fields.

    • You can change the hyperparameter value by using dot notation after you create the model. For example, Mdl.Distribution.DoF = 3.

    • To facilitate posterior sampling, bssm represents multivariate Student's t-distributed variables as an inverse-gamma scale mixture. For more details, see Latent Variance Variables of t-Distributed Errors.

For more details on distribution hyperparameters, see Distribution Hyperparameters.

Example: struct("Name","t","DoF",6) specifies a t distribution with 6 degrees of freedom for the state disturbance process.

Since R2022b

Distribution of the observation innovation process, specified as a distribution name or structure array in this table. bssm stores the value as a structure array.

DistributionNameVariable SupportStructure ArrayHyperparameter Estimation Support
Standard Gaussian"Gaussian"Multivariatestruct("Name","Gaussian")Not applicable
Student’s t"t"Multivariatestruct("Name","t","DoF",dof)Yes
Finite Gaussian mixture"mixture"Univariatestruct("Name","mixture","Weight",weight,"Mean",mu,"Variance",sigma2)No
Laplace"Laplace"Univariatestruct("Name","Laplace")Not applicable
Skew normal"skewnormal"Univariatestruct("Name","skewnormal","Delta",delta)Yes

  • The specified distribution applies to all observation innovation process variables. However, you can specify different Gaussian mixture regime means and variances among the variables.

  • For the Student's t distribution:

    • If you supply a structure array, you must specify both the "Name" and "DoF" fields.

    • You can change the hyperparameter value by using dot notation after you create the model. For example, Mdl.Distribution.DoF = 3.

    • To facilitate posterior sampling, bssm represents multivariate Student's t-distributed variables as an inverse-gamma scale mixture. For more details, see Latent Variance Variables of t-Distributed Errors.

  • For the finite Gaussian mixture distribution, bssm sets the number of regimes (mixture components) to the number of elements of weight. Therefore, to specify the nontrivial case, in which the distribution has more than one regime, you must specify at least weight.

For more details on distribution hyperparameters, see Distribution Hyperparameters.

Example: struct("Name","t","DoF",6) specifies a t distribution with 6 degrees of freedom for the state disturbance process.

Example: struct("Name","mixture","Mean",[-1 1],"Weight",[0.6 0.4]) specifies a two-regime Gaussian mixture model for the univariate observation innovation variable. The regime weights are 0.6 and 0.4. The regime mean vector is [-1 1]. The default variance for all regimes is 1.

Object Functions

estimateEstimate posterior distribution of Bayesian state-space model parameters
simulateSimulate posterior draws of Bayesian state-space model parameters
tuneTune Bayesian state-space model posterior sampler

Examples

collapse all

Create a Bayesian state-space model containing two independent, stationary, autoregressive states. The observations are the deterministic sum of the two states (in other words, the model does not impose observation errors εt). Symbolically, the system of equations is

[xt,1xt,2]=[ϕ100ϕ2][xt-1,1xt-1,2]+[σ100σ2][ut,1ut,2]

yt=[11][xt,1xt,2].

Arbitrarily assume that the prior distribution of ϕ1, ϕ2, σ12, and σ22 are independent Gaussian random variables 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 all of the following quantities:

  • A = [ϕ100ϕ2].

  • B = [σ100σ2].

  • C = [11].

  • D = 0.

  • Mean0 and Cov0 are empty arrays [], which specify the defaults.

  • StateType = [00], indicating that each state is stationary.

The paramDistribution function accepts the same vector of unknown parameters as does paramMap, but it returns the log prior density of the parameters at their current values. Specify that parameter values outside the parameter space have log prior density of -Inf.

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

Mdl = bssm(@paramMap,@priorDistribution)
Mdl = 
Mapping that defines a state-space model:
    @paramMap

Log density of parameter prior distribution:
    @priorDistribution

Mdl is a bssm 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 estimation.

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;
    Mean0 = [];         % MATLAB uses default initial state mean
    Cov0 = [];          % MATLAB uses 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

Create a time-varying, Bayesian state-space model with these characteristics:

  • From periods 1 through 250, the state equation includes stationary AR(2) and MA(1) models, respectively, and the observation model is the weighted sum of the two states.

  • From periods 251 through 500, the state model includes only the first AR(2) model.

  • μ0=[0.50.500] and Σ0 is the identity matrix.

  • The prior density is flat.

Symbolically, the state-space model is

[x1tx2tx3tx4t]=[ϕ1ϕ2001000000θ0000][x1,t-1x2,t-1x3,t-1x4,t-1]+[σ10000101][u1tu2t]yt=c1(x1t+x3t)+σ2εt.fort=1,...,250,[x1tx2t]=[ϕ1ϕ2001000][x1,t-1x2,t-1x3,t-1x4,t-1]+[σ10]u1tyt=c2x1t+σ3εt.fort=251,[x1tx2t]=[ϕ1ϕ210][x1,t-1x2,t-1]+[σ10]u1tyt=c2x1t+σ3εt.fort=252,...,500.

Write a function that specifies how the parameters theta map to the state-space model matrices, the initial state moments, and the state types. Save this code as a file named timeVariantParamMap.m on your MATLAB® path. Alternatively, open the example to access the function.

function [A,B,C,D,Mean0,Cov0,StateType] = timeVariantParamMapBayes(theta,T)
% Time-variant, Bayesian state-space model parameter mapping function
% example. This function maps the vector params to the state-space matrices
% (A, B, C, and D), the initial state value and the initial state variance
% (Mean0 and Cov0), and the type of state (StateType). From periods 1
% through T/2, the state model is a stationary AR(2) and an MA(1) model,
% and the observation model is the weighted sum of the two states. From
% periods T/2 + 1 through T, the state model is the AR(2) model only. The
% log prior distribution enforces parameter constraints (see
% flatPriorBSSM.m).
    T1 = floor(T/2);
    T2 = T - T1 - 1;
    A1 = {[theta(1) theta(2) 0 0; 1 0 0 0; 0 0 0 theta(4); 0 0 0 0]};
    B1 = {[theta(3) 0; 0 0; 0 1; 0 1]}; 
    C1 = {theta(5)*[1 0 1 0]};
    D1 = {theta(6)};
    Mean0 = [0.5 0.5 0 0];
    Cov0 = eye(4);
    StateType = [0 0 0 0];
    A2 = {[theta(1) theta(2) 0 0; 1 0 0 0]};
    B2 = {[theta(3); 0]};
    A3 = {[theta(1) theta(2); 1 0]};
    B3 = {[theta(3); 0]}; 
    C3 = {theta(7)*[1 0]};
    D3 = {theta(8)};
    A = [repmat(A1,T1,1); A2; repmat(A3,T2,1)];
    B = [repmat(B1,T1,1); B2; repmat(B3,T2,1)];
    C = [repmat(C1,T1,1); repmat(C3,T2+1,1)];
    D = [repmat(D1,T1,1); repmat(D3,T2+1,1)];
end

Write a function that specifies a joint flat prior and parameter constraints. Save this code as a file named flatPriorBSSM.m on your MATLAB path. Alternatively, open the example to access the function.

function logprior = flatPriorBSSM(theta)
% flatPriorBSSM computes the log of the flat prior density for the eight
% variables in theta (see timeVariantParamMapBayes.m). Log probabilities
% for parameters outside the parameter space are -Inf.

    % theta(1) and theta(2) are lag 1 and lag 2 terms in a stationary AR(2)
    % model. The eigenvalues of the AR(1) representation need to be within
    % the unit circle.
    evalsAR2 = eig([theta(1) theta(2); 1 0]);
    evalsOutUC = sum(abs(evalsAR2) >= 1) > 0;

    % Standard deviations of disturbances and errors (theta(3), theta(6),
    % and theta(8)) need to be positive.
    nonnegsig1 = theta(3) <= 0;
    nonnegsig2 = theta(6) <= 0;
    nonnegsig3 = theta(8) <= 0;

    paramconstraints = [evalsOutUC nonnegsig1 ...
        nonnegsig2 nonnegsig3];
    if sum(paramconstraints) > 0
        logprior = -Inf;
    else
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

Create a bssm object representing the Bayesian state-space model object. Supply the parameter-to-matrix mapping function timeVariantParamMapBayes as a function handle of solely theta by setting the time series length to 500.

numObs = 500;
Mdl = bssm(@(theta)timeVariantParamMapBayes(theta,numObs),@flatPriorBSSM)
Mdl = 
Mapping that defines a state-space model:
    @(theta)timeVariantParamMapBayes(theta,numObs)

Log density of parameter prior distribution:
    @flatPriorBSSM

This example shows how to specify that the state disturbances of a Bayesian state-space model are Student's t distributed in order to model excess kurtosis in the state equation. The example shows how to prepare the degrees of freedom parameter for posterior estimation and the example fully specifies the distribution.

Consider the Bayesian state-space model in Create Time-Invariant Bayesian State-Space Model with Known and Unknown Parameters, but assume that the state disturbances are distributed as a Student's t random variable.

Create the model by passing function handles to the local functions that represent the state-space model structure and prior distribution of the model parameters Θ. Specify that the distribution of the state disturbances is Student's t.

Mdl = bssm(@paramMap,@priorDistribution,StateDistribution="t");
Mdl.StateDistribution
ans = struct with fields:
    Name: "t"
     DoF: NaN

Mdl is a bssm model. The property Mdl.StateDistribution is a structure array specifying the distribution of the state disturbances. The field DoF is NaN by default, which means that the t-distribution degrees of freedom νu is configured for estimation with all unknown state-space model parameters Θ.

You can specify a fixed value for the degrees of freedom two ways:

  • By setting the DoF field to the positive scalar using dot notation

  • By recreating the bssm model and supplying a structure array specifying the distribution name and degrees of freedom

Specify that the degrees of freedom is 6 by using both methods.

Mdl.StateDistribution.DoF = 6
Mdl = 
Mapping that defines a state-space model:
    @paramMap

Log density of parameter prior distribution:
    @priorDistribution

Degree of freedom of t distribution in the state equation:
     6

statedist = struct("Name","t","DoF",6);
Mdl2 = bssm(@paramMap,@priorDistribution,StateDistribution=statedist)
Mdl2 = 
Mapping that defines a state-space model:
    @paramMap

Log density of parameter prior distribution:
    @priorDistribution

Degree of freedom of t distribution in the state equation:
     6

Mdl and Mdl2 are equal.

Local Functions

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;
    Mean0 = [];         % MATLAB uses default initial state mean
    Cov0 = [];          % MATLAB uses 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

To model volatility clustering, you can specify logχ12 distributed observation innovations by setting appropriate mixture weights, and regime means and variances of a 7-regime Gaussian mixture distribution.

Consider the Bayesian state-space model in Create Time-Invariant Bayesian State-Space Model with Known and Unknown Parameters, but assume that the observation-innovations process is distributed as a logχ12 random variable.

Create a structure array with the following fields and values.

  • Field Name with value "mixture"

  • Field Weight with value [0.0089 0.0541 0.1338 0.2761 0.2923 0.1494 0.0854]

  • Field Mean with value [-9.3202 -5.3145 -3.4147 -1.7097 -0.4531 0.3975 1.1925]

  • Field Variance with value [3.2793 2.4574 1.8874 1.3121 0.8843 0.5898 0.4995].^2

weight = [0.0089 0.0541 0.1338 0.2761 0.2923 0.1494 0.0854];
mu = [-9.3202 -5.3145 -3.4147 -1.7097 -0.4531 0.3975 1.1925];
sigma2 = [3.2793 2.4574 1.8874 1.3121 0.8843 0.5898 0.4995].^2;

ObsInnovDist = struct("Name","mixture","Weight",weight, ...
    "Mean",mu,"Variance",sigma2);

Create the model by passing function handles to the local functions that represent the state-space model structure and prior distribution of the model parameters Θ. Use the structure array ObsInnovDist to specify that the distribution of the observation innovations is finite Gaussian mixture with hyperparameters that define a logχ12 distribution.

Mdl = bssm(@paramMap,@priorDistribution,ObservationDistribution=ObsInnovDist);
Mdl.ObservationDistribution
ans = struct with fields:
        Name: "mixture"
      Weight: [0.0089 0.0541 0.1338 0.2761 0.2923 0.1494 0.0854]
        Mean: [-9.3202 -5.3145 -3.4147 -1.7097 -0.4531 0.3975 1.1925]
    Variance: [10.7538 6.0388 3.5623 1.7216 0.7820 0.3479 0.2495]

Mdl is a bssm model. The property Mdl.ObservationDistribution is a structure array specifying the distribution of the observation-innovations process. All distribution hyperparameters are fully specified.

Plot the distribution of the observation innovations. Compare the Gaussian mixture representation of the logχ12 and the true logχ12 distributions.

r = numel(weight);
LogChi2GMMdl = gmdistribution(mu',reshape(sigma2,1,1,r),weight);
gmPDF = @(x)arrayfun(@(x0)pdf(LogChi2GMMdl,x0),x);
logchi2PDF = @(x)((1/sqrt(2*pi))*exp((x-exp(x))/2));
figure
fplot(gmPDF,[-10,5])
hold on
fplot(logchi2PDF,"--r")
title("Log Chi-Squared Distribution: Gaussian Mixture Versus True")
legend("Gaussian Mixture","True",Location="best")

The distributions appear nearly identical.

Local Functions

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;
    Mean0 = [];         % MATLAB uses default initial state mean
    Cov0 = [];          % MATLAB uses 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

Consider a regression of the US unemployment rate onto and real gross national product (RGNP) rate, and suppose the resulting innovations are an ARMA(1,1) process. The state-space form of the relationship is

[x1,tx2,t]=[ϕθ00][x1,t-1x2,t-1]+[σ1]utyt-βZt=x1,t,

where:

  • x1,t is the ARMA process.

  • x2,t is a dummy state for the MA(1) effect.

  • yt is the observed unemployment rate deflated by a constant and the RGNP rate (Zt).

  • ut is an iid Gaussian series with mean 0 and standard deviation 1.

Load the Nelson-Plosser data set, which contains a table DataTable that has the unemployment rate and RGNP series, among other series.

load Data_NelsonPlosser

Create a variable in DataTable that represents the returns of the raw RGNP series. Because price-to-returns conversion reduces the sample size by one, prepad the series with NaN.

DataTable.RGNPRate = [NaN; price2ret(DataTable.GNPR)];
T = height(DataTable);

Create variables for the regression. Represent the unemployment rate as the observation series and the constant and RGNP rate series as the deflation data Zt.

Z = [ones(T,1) DataTable.RGNPRate];
y = DataTable.UR;

Write a function that specifies how the parameters theta map to the state-space model matrices, defers the initial state moments to the default, specifies the state types, and specifies the regression. Save this code as a file named armaDeflateYBayes.m on your MATLAB® path. Alternatively, open this example to access the function.

function [A,B,C,D,Mean0,Cov0,StateType,DeflatedY] = armaDeflateYBayes(theta,y,Z)
% Time-invariant, Bayesian state-space model parameter mapping function
% example. This function maps the vector parameters to the state-space
% matrices (A, B, C, and D), the default initial state value and the
% default initial state variance (Mean0 and Cov0), the type of state
% (StateType), and the deflated observations (DeflatedY). The log prior
% distribution enforces parameter constraints (see flatPriorDeflateY.m).
    A = [theta(1) theta(2); 0 0];
    B = [1; 1]; 
    C = [theta(3) 0];
    D = 0;
    Mean0 = [];
    Cov0 = [];
    StateType = [0 0];
    DeflatedY = y - Z*[theta(4); theta(5)];
end

Write a function that specifies a joint flat prior and parameter constraints. Save this code as a file named flatPriorDeflateY.m on your MATLAB path. Alternatively, open this example to access the function.

% Copyright 2022 The MathWorks, Inc.

function logprior = flatPriorDeflateY(theta)
% flatPriorDeflateY computes the log of the flat prior density for the five
% variables in theta (see armaDeflateYBayes.m). Log probabilities
% for parameters outside the parameter space are -Inf.

    % theta(1) and theta(2) are the AR and MA terms in a stationary
    % ARMA(1,1) model. The AR term must be within the unit circle.
    AROutUC = abs(theta(1)) >= 1;

    % The standard deviation of innovations (theta(3)) must be positive.
    nonnegsig1 = theta(3) <= 0;

    paramconstraints = [AROutUC nonnegsig1];
    if sum(paramconstraints) > 0
        logprior = -Inf;
    else
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

Create a bssm object representing the Bayesian state-space model. Specify the parameter-to-matrix mapping function as a handle to a function solely of the parameters theta.

Mdl = bssm(@(theta)armaDeflateYBayes(theta,y,Z),@flatPriorDeflateY)
Mdl = 
Mapping that defines a state-space model:
    @(theta)armaDeflateYBayes(theta,y,Z)

Log density of parameter prior distribution:
    @flatPriorDeflateY

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.

  • To specify a log χ21 distribution for the observation innovation process εt, set the ObservationDistribution to the structure array struct("Weight",weight,"Mean",mu,"Variance",sigma2), where:

    • weight is [0.0089 0.0541 0.1338 0.2761 0.2923 0.1494 0.0854].

    • mu is [-9.3202 -5.3145 -3.4147 -1.7097 -0.4531 0.3975 1.1925].

    • sigma2 is [3.2793 2.4574 1.8874 1.3121 0.8843 0.5898 0.4995].^2.

Algorithms

expand all

References

[1] Hastings, Wilfred K. "Monte Carlo Sampling Methods Using Markov Chains and Their Applications." Biometrika 57 (April 1970): 97–109. https://doi.org/10.1093/biomet/57.1.97.

[2] Metropolis, Nicholas, Rosenbluth, Arianna. W., Rosenbluth, Marshall. N., Teller, Augusta. H., and Teller, Edward. "Equation of State Calculations by Fast Computing Machines." The Journal of Chemical Physics 21 (June 1953): 1087–92. https://doi.org/10.1063/1.1699114.

Version History

Introduced in R2022a

expand all