estimate

Estimate posterior distribution of Bayesian state-space model parameters

Since R2022a

Syntax

PosteriorMdl = estimate(PriorMdl,Y,params0)

PosteriorMdl = estimate(PriorMdl,Y,params0,Name=Value)

[PosteriorMdl,estParams,EstParamCov,Summary] = estimate(___)

Description

PosteriorMdl = estimate(PriorMdl,Y,params0) returns the posterior Bayesian state-space model PosteriorMdl from combining the Bayesian state-space model prior distribution and likelihood PriorMdl with the response data Y. The input argument params0 is the vector of initial values for the unknown state-space model parameters Θ in PriorMdl.

PosteriorMdl = estimate(PriorMdl,Y,params0,Name=Value) specifies additional options using one or more name-value arguments. For example, estimate(Mdl,Y,params0,NumDraws=1e6,Thin=3,DoF=10) uses the multivariate t₁₀ distribution for the Metropolis-Hastings [1][2] proposal, draws 3e6 random vectors of parameters, and thins the sample to reduce serial correlation by discarding every 2 draws until it retains 1e6 draws.

[PosteriorMdl,estParams,EstParamCov,Summary] = estimate(___) additionally returns the following quantities using any of the input-argument combinations in the previous syntaxes:

• estParams — A vector containing the estimated parameters Θ.

• EstParamCov — The estimated variance-covariance matrix of the estimated parameters Θ.

• Summary — Estimation summary of the posterior distribution of parameters Θ. If the distribution of the state disturbances or observation innovations is non-Gaussian, Summary also includes the estimation summary of the final state values, and any estimated disturbance or innovation distribution hyperparameters.

Examples

Estimate Posterior Distribution of Time-Invariant Model

Simulate observed responses from a known state-space model, then treat the model as Bayesian and estimate the posterior distribution of the parameters by treating the state-space model parameters as unknown.

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

$$\left[\begin{array}{c}{x}_{t,1}\\ x_{t,2}\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& -0.75\end{array}\right]\left[\begin{array}{c}{x}_{t-1,1}\\ x_{t-1,2}\end{array}\right]+\left[\begin{array}{cc}1& 0\\ 0& 0.5\end{array}\right]\left[\begin{array}{c}{u}_{t,1}\\ u_{t,2}\end{array}\right]$$

$$y_t=\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{x}_{t,1}\\ x_{t,2}\end{array}\right].$$

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.

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

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

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

$$y_t=\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{x}_{t,1}\\ x_{t,2}\end{array}\right].$$

Consider a Bayesian state-space model representing the model with unknown parameters. Arbitrarily assume that the prior distribution of ${\varphi }_1$, ${\varphi }_2$, ${\sigma }_1^2$, and ${\sigma }_2^2$ 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 the following quantities:

• A = $\left[\begin{array}{cc}{\varphi }_1& 0\\ 0& {\varphi }_2\end{array}\right]$.

• B = $\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& {\sigma }_2\end{array}\right]$.

• C = $\left[\begin{array}{cc}1& 1\end{array}\right]$.

• D = 0.

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

• StateType = $\left[\begin{array}{cc}0& 0\end{array}\right]$, 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.

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 default options of estimate. Specify a random set of positive values in [0,1] to initialize the MCMC algorithm.

numParams = 4;
theta0 = rand(numParams,1);
PosteriorMdl = estimate(PriorMdl,y,theta0);

Local minimum found.

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

         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.4491  0.0905     0.3031      0.6164  
 c(2) | -0.8577  0.0606    -0.9400     -0.7365  
 c(3) |  0.9589  0.0695     0.8458      1.0699  
 c(4) |  0.4316  0.0893     0.3045      0.6023  
Proposal acceptance rate = 40.10%

PosteriorMdl.ParamMap

ans = function_handle with value:
    @paramMap

ThetaPostDraws = PosteriorMdl.ParamDistribution;
[numParams,numDraws] = size(ThetaPostDraws)

numParams = 
4

numDraws = 
1000

estimate finds an optimal proposal distribution for the Metropolis-Hastings sampler by using the tune function. Therefore, by default, estimate prints convergence information from tune. Also, estimate displays a summary of the posterior distribution of the parameters. The true values of the parameters are close to their corresponding posterior means; all are within their corresponding 95% credible intervals.

PosteriorMdl is a bssm object representing the posterior distribution.

• PosteriorMdl.ParamMap is the function handle to the function representing the state-space model structure; it is the same function as PrioirMdl.ParamMap.

• ThetaPostDraws is a 4-by-1000 matrix of draws from the posterior distribution. By default, estimate treats the first 100 draws as a burn-in sample and removes them from the matrix.

To diagnose the Markov chain induced by the Metropolis-Hastings sampler, create trace plots of the posterior parameter draws.

paramNames = ["\phi_1" "\phi_2" "\sigma_1" "\sigma_2"];
figure
h = tiledlayout(4,1);
for j = 1:numParams
    nexttile
    plot(ThetaPostDraws(j,:))
    hold on
    yline(trueTheta(j))
    ylabel(paramNames(j))
end
title(h,"Posterior Trace Plots")

The sampler eventually settles at near the true values of the parameters. In this case, the sample shows serial correlation and transient behavior. You can remedy serial correlation in the sample by adjusting the Thin name-value argument, and you can remedy transient effects by increasing the burn-in period using the BurnIn name-value argument.

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

Improve Markov Chain Convergence

Consider the model in the example Estimate Posterior Distribution of Time-Invariant Model. Improve the Markov chain convergence by adjusting sampler options.

Create a standard state-space model object ssm that represents the DGP, and then simulate a response path.

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);

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

Create a Bayesian state-space model template for estimation by passing function handles directly to paramMap and paramDistribution to bssm (the functions are in Local Functions).

Mdl = bssm(@paramMap,@priorDistribution);

Estimate the posterior distribution. Specify the simulated response path as observed responses, specify a random set of positive values in [0,1] to initialize the MCMC algorithm, and shut off all optimization displays. The plots in Estimate Posterior Distribution of Time-Invariant Model suggest that the Markov chain settles after 500 draws. Therefore, specify a burn-in period of 500 (BurnIn=500). Specify thinning the sample by keeping the first draw of each set of 30 successive draws (Thin=30). Retain 2000 random parameter vectors (NumDraws=2000).

numParams = 4;
theta0 = rand(numParams,1);
options = optimoptions("fminunc",Display="off");
PosteriorMdl = estimate(Mdl,y,theta0,Options=options, ...
    NumDraws=2000,BurnIn=500,Thin=30);

         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%

ThetaPostDraws = PosteriorMdl.ParamDistribution;

Plot trace plots and correlograms of the parameters.

paramNames = ["\phi_1" "\phi_2" "\sigma_1" "\sigma_2"];
figure
h = tiledlayout(4,1);
for j = 1:numParams
    nexttile
    plot(ThetaPostDraws(j,:))
    hold on
    yline(trueTheta(j))
    ylabel(paramNames(j))
end
title(h,"Posterior Trace Plots")

figure
h = tiledlayout(4,1);
for j = 1:numParams
    nexttile
    autocorr(ThetaPostDraws(j,:));
    ylabel(paramNames(j));
    title([]);
end
title(h,"Posterior Sample Correlograms")

The sampler quickly settles near the true values of the parameters. The sample shows little serial correlation and no transient behavior.

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

Estimate Degrees of Freedom of t-Distributed State Disturbances

Simulate observed responses from a DGP, then treat the model as Bayesian and estimate the posterior distribution of the model parameters $\Theta$ and the degrees of freedom ${\nu }_{\mathit{u}}$ of multivariate $\mathit{t}$-distributed state disturbances.

Consider the following DGP.

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

$$y_t=\left[\begin{array}{cc}1& 3\end{array}\right]\left[\begin{array}{c}{x}_{t,1}\\ x_{t,2}\end{array}\right].$$

• The true value of the state-space parameter set $\Theta =\left\{{\varphi }_1,{\varphi }_2,{\sigma }_1,{\sigma }_2\right\}=\left\{0.5,-0.75,1,0.5\right\}$.

• The state disturbances ${\mathit{u}}_{1,\mathit{t}}$ and ${\mathit{u}}_{2,\mathit{t}}$ are jointly a multivariate Student's $\mathit{t}$ random series with ${\nu }_{\mathit{u}}=5$ degrees of freedom.

Create a vector autoregression (VAR) model that represents the state equation of the DGP.

trueTheta = [0.5; -0.75; 1; 0.5];
trueDoF = 5;
phi = [trueTheta(1) 0; 0 trueTheta(2)];
Sigma = [trueTheta(3)^2 0; 0 trueTheta(4)^2];
DGP = varm(AR={phi},Covariance=Sigma,Constant=[0; 0]);

Filter a random 2-D multivariate central $\mathit{t}$ series of length 500 through the VAR model to obtain state values. Set the degrees of freedom to 5.

rng(10)  % For reproducibility
T = 500;
trueU = mvtrnd(eye(DGP.NumSeries),trueDoF,T);
X = filter(DGP,trueU);

Obtain a series of observations from the DGP by the linear combination ${\mathit{y}}_{\mathit{t}}={\mathit{x}}_{1,\mathit{e}}+3{\mathit{x}}_{2,\mathit{t}}$.

C = [1 3];
y = X*C';

Consider a Bayesian state-space model representing the model with parameters $\Theta$ and ${\nu }_{\mathit{u}}$ treated as unknown. Arbitrarily assume that the prior distribution of the parameters in $\Theta$ are independent Gaussian random variables with mean 0.5 and variance 1. Assume that the prior on the degrees of freedom ${\nu }_{\mathit{u}}$ is flat. The functions in Local Functions specify the state-space structure and prior distributions.

Create the Bayesian state-space model by passing function handles to the paramMap and priorDistribution functions to bssm. Specify that the state disturbance distribution is multivariate Student's $\mathit{t}$ with unknown degrees of freedom.

PriorMdl = bssm(@paramMap,@priorDistribution,StateDistribution="t");

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

Estimate the posterior distribution by using estimate. Specify a random set of positive values in [0,1] to initialize the MCMC algorithm. Set the burn-in period of the MCMC algorithm to 1000 draws, thin the entire MCMC sample by retaining every third draw, and set the proposal scale matrix proportionality constant to 0.25 to increase the proposal acceptance rate.

numParamsTheta = 4;
theta0 = rand(numParamsTheta,1);
PosteriorMdl = estimate(PriorMdl,y,theta0,Thin=3,BurnIn=1000,Proportion=0.25);

Local minimum found.

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

         Optimization and Tuning        
      | Params0  Optimized  ProposalStd 
----------------------------------------
 c(1) |  0.9219    0.3622      0.1151   
 c(2) |  0.9475   -0.7530      0.0454   
 c(3) |  0.2299    1.3465      0.1917   
 c(4) |  0.6759    0.5891      0.0545   
 
               Posterior Distributions              
          |   Mean     Std   Quantile05  Quantile95 
----------------------------------------------------
 c(1)     |  0.4516  0.1036     0.2641      0.5972  
 c(2)     | -0.7459  0.0376    -0.8085     -0.6863  
 c(3)     |  0.9816  0.1526     0.7430      1.2494  
 c(4)     |  0.4843  0.0402     0.4234      0.5520  
 x(1)     | -1.1901  0.9271    -2.7359      0.2539  
 x(2)     |  0.2133  0.3090    -0.2680      0.7286  
 StateDoF |  5.3980  0.7931     4.1574      6.7111  
Proposal acceptance rate = 51.10%

ThetaPostDraws = PosteriorMdl.ParamDistribution;
[numParams,numDraws] = size(ThetaPostDraws)

numParams = 
4

numDraws = 
1000

estimate displays a summary of the posterior distribution of the parameters $\Theta$ (c(1) through c(4)), the final values of the two states (x(1) and x(2)), and ${\nu }_{\mathit{u}}$ (StateDoF). The true values of the parameters are close to their corresponding posterior means; all are within their corresponding 95% credible intervals.

PosteriorMdl is a bssm object representing the posterior distribution.

• PosteriorMdl.ParamMap is the function handle to the function representing the state-space model structure. It is the same function as PrioirMdl.ParamMap.

• ThetaPostDraws is a 4-by-1000 matrix of draws from the posterior distribution of $\Theta |\mathit{y}$.

To diagnose the Markov chain induced by the Metropolis-Hastings sampler, create trace plots of the posterior parameter draws of $\Theta$.

paramNames = ["\phi_1" "\phi_2" "\sigma_1" "\sigma_2"];
figure
h = tiledlayout(4,1);
for j = 1:numParams
    nexttile
    plot(ThetaPostDraws(j,:))
    hold on
    yline(trueTheta(j))
    ylabel(paramNames(j))
end
title(h,"Posterior Trace Plots")

The sample shows some serial correlation. You can increase Thin to remedy this behavior.

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 3];
    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

Estimate Posterior of Time-Varying Model

Consider the following time-varying, state-space model for a DGP:

• 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.

• ${\mu }_0=\left[\begin{array}{cccc}0.5& 0.5& 0& 0\end{array}\right]$ and ${\Sigma }_0$ is the identity matrix.

Symbolically, the DGP is

$\begin{array}{l}\begin{array}{c}\left[\begin{array}{c}{x}_{1t}\\ x_{2t}\\ x_{3t}\\ x_{4t}\end{array}\right]=\left[\begin{array}{cccc}{\varphi }_1& {\varphi }_2& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& \theta \\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\\ x_{3,t-1}\\ x_{4,t-1}\end{array}\right]+\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& 0\\ 0& 1\\ 0& 1\end{array}\right]\left[\begin{array}{c}{u}_{1t}\\ u_{2t}\end{array}\right]\\ y_t=c_1\left(x_{1t}+x_{3t}\right)+{\sigma }_2{\epsilon }_t\end{array}fort=1,...,250,\\ \begin{array}{c}\left[\begin{array}{c}{x}_{1t}\\ x_{2t}\end{array}\right]=\left[\begin{array}{cccc}{\varphi }_1& {\varphi }_2& 0& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\\ x_{3,t-1}\\ x_{4,t-1}\end{array}\right]+\left[\begin{array}{c}{\sigma }_1\\ 0\end{array}\right]u_{1t}\\ y_t=c_2{x}_{1t}+{\sigma }_3{\epsilon }_t\end{array}fort=251,\\ \begin{array}{c}\left[\begin{array}{c}{x}_{1t}\\ x_{2t}\end{array}\right]=\left[\begin{array}{cc}{\varphi }_1& {\varphi }_2\\ 1& 0\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\end{array}\right]+\left[\begin{array}{c}{\sigma }_1\\ 0\end{array}\right]u_{1t}\\ y_t=c_2{x}_{1t}+{\sigma }_3{\epsilon }_t\end{array}fort=252,...,500,\end{array}$

where:

• The AR(2) parameters $\left\{{\varphi }_1,{\varphi }_2\right\}=\left\{0.5,-0.2\right\}$ and ${\sigma }_1=0.4$.

• The MA(1) parameter $\theta =0.3$.

• The observation equation parameters $\left\{{\mathit{c}}_1,{\mathit{c}}_2\right\}=\left\{2,3\right\}$ and $\left\{{\sigma }_2,{\sigma }_3\right\}=\left\{0.1,0.2\right\}$.

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

type timeVariantParamMapBayes.m

% Copyright 2022 The MathWorks, Inc.

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

Simulate a response path of length 500 from the model.

params = [0.5; -0.2; 0.4; 0.3; 2; 0.1; 3; 0.2];
numParams = numel(params);
numObs = 500;
[A,B,C,D,mean0,Cov0,stateType] = timeVariantParamMapBayes(params,numObs);
DGP = ssm(A,B,C,D,Mean0=mean0,Cov0=Cov0,StateType=stateType);

rng(1)  % For reproducibility
y = simulate(DGP,numObs);
plot(y)
ylabel("y")

Write a function that specifies a flat prior distribution on the state-space model parameters theta. The function returns the scalar log prior for an input set of parameters. Save this code as a file named flatPriorBSSM.m on your MATLAB® path. Alternatively, open the example to access the function.

type flatPriorBSSM.m

% Copyright 2022 The MathWorks, Inc.

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 time-varying, Bayesian state-space model that uses the structure of the DGP.

Mdl = bssm(@(params)timeVariantParamMapBayes(params,numObs),@flatPriorBSSM);

Estimate the posterior distribution. Initialize the parameter values for the MCMC algorithm with a random set of positive values in [0,0.5]. Set the proposal distribution to multivariate ${\mathit{t}}_{25}$. Set the proportionality constant of the proposal distribution scale matrix to 0.005. Shut off all displays. Return the posterior distribution, estimated posterior means of the parameters, estimated covariance matrix of the estimated parameters, and an estimation summary.

params0 = 0.5*rand(numParams,1);
options = optimoptions("fminunc",Display="off");
[PostParams,estParams,EstParamCov,Summary] = estimate(Mdl,y,params0, ...
    DoF=25,Proportion=0.005,Display=false,Options=options);
[params estParams]

ans = 8×2

    0.5000    0.5065
   -0.2000   -0.2376
    0.4000    1.0243
    0.3000    0.8592
    2.0000    0.9569
    0.1000    1.6787
    3.0000    1.2119
    0.2000    0.1781

EstParamCov

EstParamCov = 8×8

    0.0013   -0.0004    0.0006   -0.0022    0.0005   -0.0000   -0.0011   -0.0016
   -0.0004    0.0014   -0.0023   -0.0011    0.0024    0.0002    0.0025    0.0004
    0.0006   -0.0023    0.0396    0.0245   -0.0260    0.0001   -0.0436    0.0070
   -0.0022   -0.0011    0.0245    0.0745   -0.0443    0.0108   -0.0256    0.0174
    0.0005    0.0024   -0.0260   -0.0443    0.0355   -0.0069    0.0275   -0.0115
   -0.0000    0.0002    0.0001    0.0108   -0.0069    0.0038    0.0001    0.0023
   -0.0011    0.0025   -0.0436   -0.0256    0.0275    0.0001    0.0485   -0.0066
   -0.0016    0.0004    0.0070    0.0174   -0.0115    0.0023   -0.0066    0.0136

Summary

Summary=8×4 table
     Mean        Std       Quantile05    Quantile95
    _______    ________    __________    __________

    0.50653    0.035976      0.46127        0.5786 
    -0.2376     0.03725     -0.29337      -0.17496 
     1.0243     0.19899      0.75605        1.3784 
     0.8592     0.27297      0.40785        1.3276 
    0.95694     0.18852       0.6326        1.2624 
     1.6787    0.061687       1.5636        1.7739 
     1.2119     0.22027      0.83207        1.5171 
    0.17813     0.11659     0.015488       0.41109 

The sampler has some trouble estimating several parameters. You can improve posterior estimation by diagnosing the properties of the Metropolis-Hasting sample.

Estimate Posterior from Model with Deflated Response

When you work with a state-space model that contains a deflated response variable, you must have data for the predictors.

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

$$\begin{array}{l}\left[\begin{array}{c}{x}_{1,t}\\ x_{2,t}\end{array}\right]=\left[\begin{array}{cc}\varphi & \theta \\ 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\end{array}\right]+\left[\begin{array}{c}\sigma \\ 1\end{array}\right]u_t\\ y_t-\beta Z_t=x_{1,t},\end{array}$$

where:

• $$x_{1,t}$$ is the ARMA process.

• $$x_{2,t}$$ is a dummy state for the MA(1) effect.

• $$y_t$$ is the observed unemployment rate deflated by a constant and the RGNP rate ($$Z_t$$).

• $$u_t$$ 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 ${\mathit{Z}}_{\mathit{t}}$.

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

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

type armaDeflateYBayes.m

% Copyright 2022 The MathWorks, Inc.

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 = [theta(3); 1]; 
    C = [1 0];
    D = 0;
    Mean0 = [];
    Cov0 = [];
    StateType = [0 0];
    DeflatedY = y - Z*[theta(4); theta(5)];
end

Write a function that specifies a flat prior distribution on the state-space model parameters theta. The function returns the scalar log prior for an input set of parameters. Save this code as a file named flatPriorDeflateY.m on your MATLAB® path. Alternatively, open the example to access the function.

type flatPriorDeflateY.m

% 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.

numParams = 5;
Mdl = bssm(@(params)armaDeflateYBayes(params,y,Z),@flatPriorDeflateY)

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

Log density of parameter prior distribution:
    @flatPriorDeflateY

Estimate the posterior distribution. Initialize the MCMC algorithm with a random set of positive values in [0,0.5]. Set the proportionality constant to 0.01. Set a burn-in period of 2000 draws, set a thinning factor of 50, and specify retaining 1000 draws. Return an estimation summary table.

rng(1)  % For reproducibility
params0 = 0.5*rand(numParams,1);
options = optimoptions("fminunc",Display="off");
[~,~,~,Summary] = estimate(Mdl,y,params0,Proportion=0.01, ...
    BurnIn=2000,NumDraws=1000,Thin=50,Options=options,Display=false);
Summary

Summary=5×4 table
     Mean        Std       Quantile05    Quantile95
    _______    ________    __________    __________

    0.83476    0.069346     0.71916       0.94427  
     1.8213     0.24575       1.392        2.2163  
     1.7568     0.33037      1.2312        2.3438  
     7.8658      3.1097      3.0791        13.164  
    -15.168        1.92     -18.011       -11.761  

Input Arguments

PriorMdl — Prior Bayesian state-space model

bssm model object

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

The function handles of the properties PriorMdl.ParamDistribution and PriorMdl.ParamMap determine the prior and the data likelihood, respectively.

Y — Observed response data

numeric matrix | cell vector of numeric vectors

Observed response data, from which estimate 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

params0 — Initial values for parameters Θ

numeric vector

Initial parameter values for the parameters Θ, specified as a numParams-by-1 numeric vector. Elements of params0 must correspond to the elements of the first input arguments of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

Data Types: double

Name-Value Arguments

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.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: estimate(Mdl,Y,params0,NumDraws=1e6,Thin=3,DoF=10) uses the multivariate t₁₀ distribution for the Metropolis-Hastings proposal, draws 3e6 random vectors of parameters, and thins the sample to reduce serial correlation by discarding every 2 draws until it retains 1e6 draws.

Kalman Filter Options

Univariate — Univariate treatment of multivariate series flag

false (default) | true

Univariate treatment of a multivariate series flag, specified as a value in this table.

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 and PriorMdl.ParamDistribution.

Example: Univariate=true

Data Types: logical

SquareRoot — Square root filter method flag

false (default) | true

Square root filter method flag, specified as a value in this table.

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

If you suspect that 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

Posterior Sampling Options

NumDraws — Number of Metropolis-Hastings sampler draws

1000 (default) | positive integer

Number of Metropolis-Hastings sampler draws from the posterior distribution Π(θ|Y), specified as a positive integer.

Example: NumDraws=1e7

Data Types: double

BurnIn — Number of draws to remove from beginning of sample

100 (default) | nonnegative scalar

Number of draws to remove from the beginning of the sample to reduce transient effects, specified as a nonnegative scalar. For details on how estimate reduces the full sample, see Algorithms.

Tip

To help you specify the appropriate burn-in period size:

1. Determine the extent of the transient behavior in the sample by setting the BurnIn name-value argument to 0.

2. Simulate a few thousand observations by using simulate.

3. Create trace plots.

Example: BurnIn=1000

Data Types: double

Thin — Adjusted sample size multiplier

1 (no thinning) (default) | positive integer

Adjusted sample size multiplier, specified as a positive integer.

The actual sample size is BurnIn + NumDraws*Thin. After discarding the burn-in, estimate discards every Thin – 1 draws, and then retains the next draw. For more details on how estimate reduces the full sample, see Algorithms.

Tip

To reduce potential large serial correlation in the posterior sample, or to reduce the memory consumption of the output sample, specify a large value for Thin.

Example: Thin=5

Data Types: double

DoF — Proposal distribution degrees of freedom

Inf (default) | positive scalar

Proposal distribution degrees of freedom for parameter updates using the Metropolis-Hastings sampler, specified as a value in this table.

Value	Metropolis-Hastings Proposal Distribution
Positive scalar	Multivariate t with DoF degrees of freedom
Inf	Multivariate normal

The following options specify other aspects of the proposal distribution:

• Proportion — Optional proportionality constant that scales Proposal

• Center — Optional expected value

• StdDoF — Optional proposal standard deviation of the degrees of freedom parameter of the t-distributed state disturbance or observation innovation process.

Example: DoF=10

Data Types: double

Proportion — Proposal scale matrix proportionality constant

1 (default) | positive scalar

Proposal scale matrix proportionality constant, specified as a positive scalar.

Tip

For higher proposal acceptance rates, experiment with relatively small values for Proportion.

Example: Proportion=1

Data Types: double

Center — Proposal distribution center

[] (default) | numeric vector

Proposal distribution center, specified as a value in this table.

Value	Description
numParams-by-1 numeric vector	estimate uses the independence Metropolis-Hastings sampler. Center is the center of the proposal distribution.
[] (empty array)	estimate uses the random-walk Metropolis-Hastings sampler. The center of the proposal density is the current state of the Markov chain.

Example: Center=ones(10,1)

Data Types: double

StdDoF — Proposal standard deviation for degrees of freedom parameter

0.1 (default) | positive numeric scalar

Proposal standard deviation for the degrees of freedom parameter of the t-distributed state disturbance or observation innovation process, specified as a positive numeric scalar.

StdDoF applies to the corresponding process when PriorMdl.StateDistribution.Name is "t" or PriorMdl.ObservationDistribution.Name is "t", and their associated degrees of freedom are estimable (the DoF field is NaN or a function handle). For example, if the following conditions hold, StdDof applies only to the t-distribution degrees of freedom of the state disturbance process:

• PriorMdl.StateDistribution.Name is "t".

• PriorMdl.StateDistribution.DoF is NaN.

• The limiting distribution of the observation innovations is Gaussian (PriorMdl.ObservationDistribution.Name is "Gaussian" or PriorMdl.ObservationDistribution is struct("Name","t","DoF",Inf).

Example: StdDoF=0.5

Data Types: double

Proposal Tuning Options

Hessian — Hessian approximation method

"difference" (default) | "diagonal" | "opg" | "optimizer" | character vector

Hessian approximation method for the Metropolis-Hastings proposal distribution scale matrix, specified as a value in this table.

Value	Description
"difference"	Finite differencing
"diagonal"	Diagonalized result of finite differencing
"opg"	Outer product of gradients, ignoring the prior distribution
"optimizer"	Posterior distribution optimized by fmincon or fminunc. Specify optimization options by using the Options name-value argument.

Tip

The Hessian="difference" setting can be computationally intensive and inaccurate, and the resulting scale matrix can be nonnegative definite. Try one of the other options for better results.

Example: Hessian="opg"

Data Types: char | string

Lower — Parameter lower bounds

[] (default) | numeric vector

Parameter lower bounds when computing the Hessian matrix (see Hessian), specified as a numParams-by-1 numeric vector.

Lower(j) specifies the lower bound of parameter theta(j), the first input argument of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

The default value [] specifies no lower bounds.

Note

Lower does not apply to posterior simulation. To apply parameter constraints on the posterior, code them in the log prior distribution function PriorMdl.ParamDistribution by setting the log prior of values outside the distribution support to -Inf.

Example: Lower=[0 -5 -1e7]

Data Types: double

Upper — Parameter upper bounds

[] (default) | numeric vector

Parameter lower bounds when computing the Hessian matrix (see Hessian), specified as a numParams-by-1 numeric vector.

Upper(j) specifies the upper bound of parameter theta(j), the first input argument of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

The default value [] specifies no upper bounds.

Note

Upper does not apply to posterior simulation. To apply parameter constraints on the posterior, code them in the log prior distribution function PriorMdl.ParamDistribution by setting the log prior of values outside the distribution support to -Inf.

Example: Upper=[5 100 1e7]

Data Types: double

Options — Optimization options

optimoptions optimization controller

Optimization options for the setting Hessian="optimizer", specified as an optimoptions optimization controller. Options replaces default optimization options of the optimizer. For details on altering default values of the optimizer, see the optimization controller optimoptions, the constrained optimization function fmincon, or the unconstrained optimization function fminunc in Optimization Toolbox™.

For example, to change the constraint tolerance to 1e-6, set options = optimoptions(@fmincon,ConstraintTolerance=1e-6,Algorithm="sqp"). Then, pass Options by using Options=options.

By default, estimate uses the default options of the optimizer.

Simplex — Simplex search flag

true (default) | false

Simplex search flag to improve initial parameter values, specified as a value in this table.

Value	Description
true	Apply simplex search method to improve initial parameter values for proposal optimization. For more details, see fminsearch Algorithm.
false	Does not apply simplex search method.

estimate applies simplex search when the numerical optimization exit flag is not positive.

Example: Simplex=false

Data Types: logical

Display — Proposal tuning results display flag

true (default) | false

Proposal tuning results display flag, specified as a value in this table.

Value	Description
true	Displays tuning results
false	Suppresses tuning results display

Example: Display=false

Data Types: logical

Output Arguments

PosteriorMdl — Posterior Bayesian state-space model

bssm model object

Posterior Bayesian state-space model, returned as a bssm model object.

The property PosteriorMdl.ParamDistribution is a numParams-by-NumDraws sample from the posterior distribution of Θ|y; estimate obtains the sample by using the Metropolis-Hastings sampler. PosteriorMdl.ParamDistribution(j,:) contains a NumDraws sample of the parameter theta(j), where theta corresponds to Θ, and it is the first input argument of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

The properties PosteriorMdl.ParaMap and PrioirMdl.ParamMap are equal.

Tip

To obtain posterior draws of distribution hyperparameters, such as the degrees of freedom parameters ν_u and ν_ε of corresponding t-distributed errors, use simulate and specify an output function, using the OutputFunction name-value argument, that returns the draws after each MCMC iteration.

estParams — Posterior means of Θ

numeric vector

Posterior means of Θ, returned as a numParams-by-1 numeric vector. estParams(j) contains the mean of the NumDraws posterior sample of the parameter theta(j).

EstParamCov — Variance-covariance matrix of estimates of Θ

numeric matrix

Variance-covariance matrix of estimates of Θ, returned as a numParams-by-numParams numeric matrix. EstParamCov(j,k) contains the estimated covariance of the parameter estimates of theta(j) and theta(k), based on the NumDraws posterior sample of those parameters.

Summary — Summary of all Bayesian estimators

structure array

Summary of all Bayesian estimators, returned as a table.

Rows of the table correspond to elements of Θ, final state values, and estimated state disturbance and observation innovation hyperparameters:

• For the first numParams rows of the table, row j contains the posterior estimation summary of theta(j).

• If the distribution of the state disturbances or observation innovations is non-Gaussian, the following statements hold:

  • Rows numParams + k contains the posterior estimation summary of the final value of state k.

  • Rows after numParams + m_T contain the posterior estimation summary of any estimated state disturbance hyperparameters, followed by any estimated observation innovation hyperparameters. m_T is the number of states in the final period T.

Columns contain the posterior mean Mean, standard deviation Std, 5% percentile Quantile05, and 95% percentile Quantile95.

Algorithms

• The Kalman filter accommodates missing data by not updating filtered state estimates corresponding to missing observations. In other words, suppose there is a missing observation at period t. Then, the state forecast for period t based on the previous t – 1 observations and filtered state for period t are equivalent.

• estimate uses the tune function to create the proposal distribution for the Metropolis-Hastings sampler. You can tune the sampler by using the proposal tuning options.

• When estimate tunes the proposal distribution, the optimizer that estimate uses to search for the posterior mode before computing the Hessian matrix depends on your specifications.

  • estimate uses fminunc by default.

  • estimate uses fmincon when you perform constrained optimization by specifying the Lower or Upper name-value argument.

  • estimate further refines the solution returned by fminunc or fmincon by using fminsearch when Simplex=true, which is the default.

• This figure shows how estimate reduces the sample by using the values of NumDraws, Thin, and BurnIn. Rectangles represent successive draws from the distribution. estimate removes the white rectangles from the sample. The remaining NumDraws black rectangles compose the sample.

  

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