# estimate

Estimate posterior distribution of Bayesian state-space model parameters

## Syntax

PosteriorMdl = estimate(PriorMdl,Y,params0)
PosteriorMdl = estimate(PriorMdl,Y,params0,Name=Value)
[PosteriorMdl,estParams,EstParamCov,Summary] = estimate(___)

## Description

example

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.

example

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

example

[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 parametersEstParamCov — The estimated variance-covariance matrix of the estimated parametersSummary — Estimation summary of the posterior distribution 

## Examples

collapse all

Simulate observed responses from a known state-space model, then treat the model as Bayesian and estimate the posterioir 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 the 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 Kalman filter.

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

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 the Bayesian state-space model by passing function handles directly to paramMap and paramDistribution to bssm (the functions are in Local Functions).

Mdl = bssm(@paramMap,@priorDistribution)
Mdl = Mapping that defines a state-space model: @paramMap Log density of the prior distribution: @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 Kalman filter, 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

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}\phantom{\rule{0.2777777777777778em}{0ex}}for\phantom{\rule{0.2777777777777778em}{0ex}}t=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}\phantom{\rule{0.2777777777777778em}{0ex}}for\phantom{\rule{0.2777777777777778em}{0ex}}t=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}\phantom{\rule{0.2777777777777778em}{0ex}}for\phantom{\rule{0.2777777777777778em}{0ex}}t=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\}$.

Simulate a response path of length 500 from the model. The function timeVariantParamMapBayes.m, stored in matlabroot/examples/econ/main, specifies the state-space model structure.

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") Create a time-varying, Bayesian state-space model that uses the structure of the DGP, but all parameters are unknown and the prior density is flat (the function flatPriorBSSM.m in matlabroot/examples/econ/main is the prior density).

Create a bssm object representing the Bayesian state-space model object.

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

Estimate the posterior distribution. Initialize the parameter values for the Kalman filter 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.

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;

The functions armaDeflateYBayes.m and flatPriorDeflateY.m in matlabroot/examples/econ/main specify the state-space model structure and likelihood, and the flat prior distribution.

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 the prior distribution: @flatPriorDeflateY 

Estimate the posterior distribution. Initialize the Kalman filter 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

collapse all

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.

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

Initial parameter values, 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 Mdl.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 t10 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

collapse all

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

ValueDescription
trueApplies the univariate treatment of a multivariate series, also known as sequential filtering
falseDoes 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, DtDt' must be diagonal, where Dt (t = 1, ..., T) is the output coefficient matrix D of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

Example: Univariate=true

Data Types: logical

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

ValueDescription
trueApplies the square root filter method for the Kalman filter
falseDoes 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

collapse all

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

Example: NumDraws=1e7

Data Types: double

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

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

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

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

ValueMetropolis-Hastings Proposal Distribution
Positive scalarMultivariate t with Dof degrees of freedom
InfMultivariate normal

The following options specify other aspects of the proposal distribution:

Example: Dof=10

Data Types: double

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

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

ValueDescription
numParams-by-1 numeric vectorestimate 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

Proposal Tuning Options

collapse all

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

ValueDescription
"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

Parameter lower bounds when computing the Hessian matrix (see Hessian), specified as a numParams-by-1 numeric vector. If you specify Lower, estimate uses

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

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

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 search flag to improve initial parameter values, specified as a value in this table.

ValueDescription
trueApply simplex search method to improve initial parameter values for proposal optimization. For more details, see fminsearch Algorithm.
falseDoes not apply simplex search method.

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

Example: Simplex=false

Data Types: logical

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

ValueDescription
trueDisplays tuning results
falseSuppresses tuning results display

Example: Display=false

Data Types: logical

## Output Arguments

collapse all

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; estimate obtains the sample by using the Metropolis-Hastings sampler. PosteriorMdl.ParamDistribution(j,:) contains a NumDraws sample of the parameter theta(j), where theta is the first input argument of PriorMdl.ParamMap and PriorMdl.ParamDistribution.

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

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

Variance-covariance matrix of the parameter estimates, 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).

Summary of Bayesian estimators, returned as a table. Row j of the table contains the posterior estimation summary of theta(j). Columns contain the posterior mean Mean, standard deviation Std, 5% percentile Quantile05, and 95% percentile Quantile95.

## Algorithms

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

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

 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.