estimate
Estimate posterior distribution of Bayesian state-space model parameters
Since R2022a
Syntax
Description
returns the posterior Bayesian state-space model PosteriorMdl
= estimate(PriorMdl
,Y
,params0
)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
.
specifies additional options using one or more name-value arguments. For example,
PosteriorMdl
= estimate(PriorMdl
,Y
,params0
,Name=Value
)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.
[
additionally returns the following quantities using any of the input-argument combinations
in the previous syntaxes:PosteriorMdl
,estParams
,EstParamCov
,Summary
]
= estimate(___)
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).
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
Consider a Bayesian state-space model representing the model with unknown parameters. Arbitrarily assume that the prior distribution of , , , and 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
= .B
= .C
= .D
= 0.Mean0
andCov0
are empty arrays[]
, which specify the defaults.StateType
= , 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 asPrioirMdl.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 and the degrees of freedom of multivariate -distributed state disturbances.
Consider the following DGP.
The true value of the state-space parameter set .
The state disturbances and are jointly a multivariate Student's random series with 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 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 .
C = [1 3]; y = X*C';
Consider a Bayesian state-space model representing the model with parameters and treated as unknown. Arbitrarily assume that the prior distribution of the parameters in are independent Gaussian random variables with mean 0.5 and variance 1. Assume that the prior on the degrees of freedom 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 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 (c(1)
through c(4)
), the final values of the two states (x(1)
and x(2)
), and (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 asPrioirMdl.ParamMap
.ThetaPostDraws
is a 4-by-1000 matrix of draws from the posterior distribution of .
To diagnose the Markov chain induced by the Metropolis-Hastings sampler, create trace plots of the posterior parameter draws of .
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.
and is the identity matrix.
Symbolically, the DGP is
where:
The AR(2) parameters and .
The MA(1) parameter .
The observation equation parameters and .
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 . 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
where:
is the ARMA process.
is a dummy state for the MA(1) effect.
is the observed unemployment rate deflated by a constant and the RGNP rate ().
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 .
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
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 ofY
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 ofPriorMdl.ParamMap
,C{t}
, andD{t}
must be consistent with the matrix inY{t}
for all periods. The last cell ofY
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 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.
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,
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
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
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:
Determine the extent of the transient behavior in the sample by setting the
BurnIn
name-value argument to0
.Simulate a few thousand observations by using
simulate
.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 scalesProposal
Center
— Optional expected valueStdDoF
— 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
isNaN
.The limiting distribution of the observation innovations is Gaussian (
PriorMdl.ObservationDistribution.Name
is"Gaussian"
orPriorMdl.ObservationDistribution
isstruct("Name","t","DoF",Inf)
.
Example: StdDoF=0.5
Data Types: double
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(
specifies the lower bound of parameter j
)theta(
, the first input argument of j
)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(
specifies the upper bound of parameter j
)theta(
, the first input argument of j
)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(
contains a j
,:)NumDraws
sample of the parameter
theta(
, where
j
)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(
contains the mean of
the j
)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(
contains the estimated covariance of the parameter estimates of
j
,k
)theta(
and
j
)theta(
, based on the
k
)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, rowj
contains the posterior estimation summary oftheta(
.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 statek
.Rows after
numParams
+ mT contain the posterior estimation summary of any estimated state disturbance hyperparameters, followed by any estimated observation innovation hyperparameters. mT 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 thetune
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 thatestimate
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 ofNumDraws
,Thin
, andBurnIn
. Rectangles represent successive draws from the distribution.estimate
removes the white rectangles from the sample. The remainingNumDraws
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
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)