Main Content

estimate

Estimate posterior distribution of Bayesian state-space model parameters

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

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 parameters

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

  • Summary — 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).

[xt,1xt,2]=[0.500-0.75][xt-1,1xt-1,2]+[1000.5][ut,1ut,2]

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

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

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

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

Consider a Bayesian state-space model representing the model with unknown parameters. 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 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.

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

Figure contains 4 axes objects. Axes object 1 contains 2 objects of type line, constantline. Axes object 2 contains 2 objects of type line, constantline. Axes object 3 contains 2 objects of type line, constantline. Axes object 4 contains 2 objects of type line, constantline.

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 contains 4 axes objects. Axes object 1 contains 2 objects of type line, constantline. Axes object 2 contains 2 objects of type line, constantline. Axes object 3 contains 2 objects of type line, constantline. Axes object 4 contains 2 objects of type line, constantline.

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

Figure contains 4 axes objects. Axes object 1 contains 4 objects of type stem, line. Axes object 2 contains 4 objects of type stem, line. Axes object 3 contains 4 objects of type stem, line. Axes object 4 contains 4 objects of type stem, line.

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.

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

Symbolically, the DGP is

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

where:

  • The AR(2) parameters {ϕ1,ϕ2}={0.5,-0.2} and σ1=0.4.

  • The MA(1) parameter θ=0.3.

  • The observation equation parameters {c1,c2}={2,3} and {σ2,σ3}={0.1,0.2}.

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

Figure contains an axes object. The axes object contains an object of type line.

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

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

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

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

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:

  • Proportion — Optional proportionality constant that scales Proposal

  • Center — Optional expected value

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.

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

    Sample reduced by for values of NumDraws, Thin, and BurnIn

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