estimate

Class: ssm

Estimate state-space model parameters

Syntax

  • EstMdl = estimate(Mdl,Y,params0) example
  • EstMdl = estimate(Mdl,Y,params0,Name,Value) example
  • [EstMdl,estParams,EstParamCov,logL,Output] = estimate(___) example

Description

example

EstMdl = estimate(Mdl,Y,params0) estimates the parameters of Mdl using the Kalman filter and maximum likelihood, where:

  • Mdl is a state-space model (ssm).

  • Y is the observed response series.

  • params0 is the vector of initial values for unknown parameters.

The software returns EstMdl, which is the estimated state-space model (ssm) that stores the estimated coefficient matrices and initial state means, covariances, and distributions.

  • For explicitly specified state-space models (that is, you specified Mdl without using a parameter-to-matrix mapping function), the software estimates all NaN values in the coefficient matrices (Mdl.A, Mdl.B, Mdl.C, and Mdl.D) and the initial state means and covariance matrix (Mdl.Mean0 and Mdl.Cov0).

  • For implicitly specified state-space models (that is, you specified Mdl using a parameter-to-matrix mapping function, which the software stores in Mdl.ParamMap), you define the model and the location of the unknown parameters using the parameter-to-matrix mapping function. Implicitly specify a state-space model to estimate complex models, impose parameter constraints, and estimate initial states. The parameter-to-mapping function can also accommodate additional output arguments.

example

EstMdl = estimate(Mdl,Y,params0,Name,Value) estimates the state-space model Mdl with additional options specified by one or more Name,Value pair arguments.

For example, pass in predictor data to include a linear regression component to the observation equation, control how the results appear in the Command Window, and indicate which estimation method to use for the parameter covariance matrix.

example

[EstMdl,estParams,EstParamCov,logL,Output] = estimate(___) additionally returns:

  • estParams, a vector containing the estimated parameters

  • EstParamCov, the estimated variance-covariance matrix of the estimated parameters

  • logL, the optimized loglikelihood value

  • Output, optimization diagnostic information structure

using any of the input arguments in the previous syntaxes.

Tips

Constrained likelihood objective function maximization

  • You can specify any combination of linear inequality, linear equality, and upper and lower bound constraints on the parameters.

  • If a parameter is unbounded below, then set 'lb',-Inf.

  • If a parameter is unbounded above, then set 'ub',Inf.

  • It is good practice to avoid equality and inequality constraints during optimization. For example, if you want to constrain the parameter w to be positive, then implicitly specify the state-space model using a parameter-to-matrix mapping function, set w = exp(s) within the function, and use unconstrained optimization to estimate s. Subsequently, s can assume any real value, but w must be positive.

Predictors and corresponding coefficients

  • The state-space model Mdl does not store the predictors (Zt) nor their corresponding regression coefficients (β). Supply the predictors and their corresponding coefficients wherever necessary using the appropriate name-value pair arguments.

  • The predictor series serve as observation deflators. Subsequently, the deflated data set is YtZtβ, where:

    • Zt=(z1tz2tzdt)., that is, Z is a T-byd matrix.

    • zjt is the period t value of predictor j.

    • β is a d-by-n matrix of regression coefficients.

  • To include an overall mean to the observation model, include a column of 1s in Zt.

  • If you want to account for predictor effects when you simulate (simulate), then you must deflate the observations manually. To deflate the observations, use Wt=YtZtβ^.

  • If the state equation requires predictors, then expand the states by the constant 1 and the predictors.

  • If the regression model is complex, then consider implicitly defining the state space model. For example, define the parameter-to-matrix mapping function using the following syntax pattern.

    function [A,B,C,D,Mean0,Cov0,StateType,DeflateY] = ParamMap(params,Y,Z)
    		...
    		DeflateY = Y - exp(params(9) + params(10)*Z);
    		...
    end

    In this example, Y is the matrix of observations and Z is the matrix of predictors. The function returns DeflateY, which is the matrix of deflated observations. Specify Y and Z in the MATLAB® Workspace before, and then pass ParamMap to ssm using the following syntax pattern.

    Mdl = ssm(@(params)ParamMap(params,Y,Z))

    This is also useful if each response series requires a distinct set of predictors.

  • If the state equation requires predictors, then any one of the following:

    • Expanding the states by including the constant 1 state.

    • Expanding the states by including the predictors.

  • If the model is time varying with respect the observed responses, then the software does not support including predictors. If the observation vectors among different periods vary in length, then the software cannot determine which coefficients to use to deflate the observed responses.

Additional Tips

  • The software accommodates missing data. Indicate missing data using NaN values in the observed responses (Y).

  • It is good practice to check the convergence status of the optimization routine by displaying Output.ExitFlag.

  • If the optimization algorithm does not converge, then you can increase the number of iterations using the 'Options' name-value pair argument.

  • If the optimization algorithm does not converge, then consider using refine, which might help you obtain better initial parameter values for optimization.

Input Arguments

expand all

Mdl — State-space modelssm model

State-space model containing unknown parameters, specified as an ssm model returned by ssm.

Mdl does not store observed responses or predictor data. Supply the data wherever necessary using the appropriate input and name-value pair arguments.

Y — Observed response datacell vector of numeric vectors | matrix

Observed response data to which Mdl is fit, specified as a cell vector of numeric vectors or a matrix.

  • If Mdl is time invariant with respect to the observation equation, then Y is a T-by-n matrix, where each row corresponds to a period and each column corresponds to a particular observation in the model. Therefore, T is the sample size and m is the number of observations per period. The last row of Y contains the latest observations.

  • If Mdl is time varying with respect to the observation equation, then Y is a T-by-1 cell vector. Each element of the cell vector corresponds to a period and contains an nt-dimensional vector of observations for that period. The corresponding dimensions of the coefficient matrices in Mdl.C{t} and Mdl.D{t} must be consistent with the matrix in Y{t} for all periods. The last cell of Y contains the latest observations.

  • Suppose you created Mdl implicitly by specifying a parameter-to-matrix mapping function, and the function has input arguments for the observed responses or predictors (see the ParamMap name-value pair argument of ssm), the function then establishes a link to observed responses and the predictor data. The link overrides the value of Y.

Data Types: double | cell

params0 — Unknown parameter initial valuesnumeric vector

Unknown parameter initial values for numerical maximum likelihood estimation, specified as a numeric vector.

The elements of params0 correspond to the unknown parameters in the state-space model matrices A, B, C, and D, and, optionally, the initial state mean Mean0 and covariance matrix Cov0.

  • If you created Mdl explicitly (that is, by specifying the matrices without a parameter-to-matrix mapping function), then the software maps the elements of params to NaNs in the state-space model matrices and initial state values. The software searches for NaNs column-wise following the order A, B, C, D, Mean0, and Cov0.

  • If you created Mdl implicitly (that is, by specifying the matrices with a parameter-to-matrix mapping function), then you must set initial parameter values for the state-space model matrices, initial state values, and state types within the parameter-to-matrices mapping function.

Data Types: double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Estimation Options

'Beta0' — Initial values of regression coefficientsmatrix

Initial values of regression coefficients, specified as the comma-separated pair consisting of 'Beta0' and a matrix.

The number of rows in Beta0 must equal the number of columns in Y. Beta0 and Predictors must have the same number of columns.

By default, Beta0 is the ordinary least-squares estimate of Y onto Predictors.

Data Types: double

'CovMethod' — Asymptotic covariance estimation method'OPG' (default) | 'Hessian' | 'Sandwich'

Asymptotic covariance estimation method, specified as the comma-separated pair consisting of 'CovMethod' and a string.

Set CovMethod using a value in this table.

ValueDescription
'Hessian'Negative, inverted Hessian matrix
'OPG'Outer product of gradients (OPG)
'Sandwich'Both Hessian and OPG

Example: 'CovMethod','Sandwich'

Data Types: char

'Display' — Command Window display option'params' (default) | 'diagnostics' | 'full' | 'iter' | 'off' | cell vector of strings

Command Window display option, specified as the comma-separated pair consisting of 'Display' and a string or cell vector of strings.

Set Display using any combination of values in this table.

Valueestimate Displays
'diagnostics'Optimization diagnostics
'full'Maximum likelihood parameter estimates, standard errors, t statistics, iterative optimization information, and optimization diagnostics
'iter'Iterative optimization information
'off'Nothing in the Command Window
'params'Maximum likelihood parameter estimates, standard errors, and t statistics

For example,

  • To run a simulation where you are fitting many models, and therefore want to suppress all output, use 'Display','off'.

  • To display all estimation results and the optimization diagnostics, use 'Display',{'params','diagnostics'}.

Data Types: char | cell

'Options' — Optimization optionsoptimoptions optimization controller

Optimization options, specified as the comma-separated pair consisting of 'Options' and 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, suppose that you want to change the constraint tolerance to 1e-6. Set Options = optimoptions(@fmincon,'TolCon',1e-6,'Algorithm','sqp') and then pass Options into estimate using 'Options',Options.

By default:

  • For constrained optimization, estimate maximizes the likelihood objective function using fmincon and its default options, but sets 'Algorithm','interior-point'.

  • For unconstrained optimization, estimate maximizes the likelihood objective function using fminunc and its default options, but sets 'Algorithm','quasi-newton'.

'Predictors' — Predictor variables in state-space model observation equation[] (default) | matrix

Predictor variables in the state-space model observation equation, specified as the comma-separated pair consisting of 'Predictors' and a matrix. The columns of Predictors correspond to individual predictor variables. Predictors must have T rows, where row t corresponds to the observed predictors at period t (Zt) in the expanded observation equation

ytZtβ=Cxt+Dut.

That is, the software deflates the observations using the regression component. β is the time-invariant vector of regression coefficients that the software estimates with all other parameters.

If there are n observations per period, then the software regresses all predictor series onto each observation. Subsequently, the software returns d-by-n matrix of fitted regression coefficient vectors for each observation series, where d is the number of columns of Predictors.

If you specify Predictors, then Mdl must be time invariant. Otherwise, the software returns an error.

By default, the software excludes a regression component from the state-space model.

Data Types: double

'SquareRoot' — Square root filter method flagfalse (default) | true

Square root filter method flag, specified as the comma-separated pair consisting of 'SquareRoot' and true or false.

Example: 'SquareRoot',true

Data Types: logical

'Tolerance' — Forecast uncertainty threshold0 (default) | nonnegative scalar

Forecast uncertainty threshold, specified as the comma-separated pair consisting of 'Tolerance' and a nonnegative scalar.

If the forecast uncertainty for a particular observation is less than Tolerance during numerical estimation, then the software removes the uncertainty corresponding to the observation from the forecast covariance matrix before its inversion.

It is best practice to set Tolerance to a small number, for example, le-15, to overcome numerical obstacles during estimation.

Example: 'Tolerance',le-15

Data Types: double

'Univariate' — Univariate treatment of multivariate series flagfalse (default) | true

Univariate treatment of a multivariate series flag, specified as the comma-separated pair consisting of 'Univariate' and true or false.

DtDt' must be diagonal, where Dt is one of the following:

  • The matrix D{t} in a time-varying state-space model

  • The matrix D in a time-invariant state-space model

Example: 'Univariate',true

Data Types: logical

Constrained Optimization Options for fmincon

'Aeq' — Linear equality constraint parameter transformermatrix

Linear equality constraint parameter transformer for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'Aeq' and a matrix.

If you specify Aeq and beq, then estimate maximizes the likelihood objective function using the equality constraint Aeqθ=beq, where θ is a vector containing every Mdl parameter.

The number of rows of Aeq is the number of constraints, and the number of columns is the number of parameters that the software estimates. Order the columns of Aeq by Mdl.A, Mdl.B, Mdl.C, Mdl.D, Mdl.Mean0, Mdl.Cov0, and the regression coefficient (if the model has one).

Specify Aeq and beq together, otherwise estimate returns an error.

Aeq directly corresponds to the input argument Aeq of fmincon, not to the state-transition coefficient matrix Mdl.A.

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

'Aineq' — Linear inequality constraint parameter transformermatrix

Linear inequality constraint parameter transformer for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'Aineq' and a matrix.

If you specify Aineq and bineq, then estimate maximizes the likelihood objective function using the inequality constraint Aineqθbineq, where θ is a vector containing every Mdl parameter.

The number of rows of Aineq is the number of constraints, and the number of columns is the number of parameters that the software estimates. Order the columns of Aineq by Mdl.A, Mdl.B, Mdl.C, Mdl.D, Mdl.Mean0, Mdl.Cov0, and the regression coefficient (if the model has one).

Specify Aineq and bineq together, otherwise estimate returns an error.

Aineq directly corresponds to the input argument A of fmincon, not to the state-transition coefficient matrix Mdl.A.

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

Data Types: double

'beq' — Linear equality constraints of transformed parametersnumeric vector

Linear equality constraints of the transformed parameters for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'beq' and a numeric vector.

If you specify Aeq and beq, then estimate maximizes the likelihood objective function using the equality constraint Aeqθ=beq, where θ is a vector containing every Mdl parameter..

Specify Aeq and beq together, otherwise estimate returns an error.

beq directly corresponds to the input argument beq of fmincon, and is not associated with any component of Mdl.

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

Data Types: double

'bineq' — Linear inequality constraint upper boundsnumeric vector

Linear inequality constraint upper bounds of the transformed parameters for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'bineq' and a numeric vector.

If you specify Aineq and bineq, then estimate maximizes the likelihood objective function using the inequality constraint Aineqθbineq, where θ is a vector containing every Mdl parameter.

Specify Aineq and bineq together, otherwise estimate returns an error.

bineq directly corresponds to the input argument b of fmincon, and is not associated with any component of Mdl.

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

Data Types: double

'lb' — Lower bounds of parametersnumeric vector

Lower bounds of the parameters for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'lb' and a numeric vector.

If you specify lb and ub, then estimate maximizes the likelihood objective function subject tolbθub, where θ is a vector containing every Mdl parameter.

Order the elements of lb by Mdl.A, Mdl.B, Mdl.C, Mdl.D, Mdl.Mean0, Mdl.Cov0, and the regression coefficient (if the model has one).

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

Data Types: double

'ub' — Upper bounds of parametersnumeric vector

Upper bounds of the parameters for constrained likelihood objective function maximization, specified as the comma-separated pair consisting of 'ub' and a numeric vector.

If you specify lb and ub, then estimate maximizes the likelihood objective function subject tolbθub, where θ is a vector every Mdl parameter.

Order the elements of ub by Mdl.A, Mdl.B, Mdl.C, Mdl.D, Mdl.Mean0, Mdl.Cov0, and the regression coefficient (if the model has one).

By default, if you did not specify any constraint (linear inequality, linear equality, or upper and lower bound), then estimate maximizes the likelihood objective function using unconstrained maximization.

Data Types: double

Output Arguments

expand all

EstMdl — State-space model containing parameter estimatesssm model

State-space model containing the parameter estimates, returned as an ssm model.

estimate uses maximum likelihood to calculate all parameter estimates. EstMdl stores the parameter estimates in the coefficient matrices (EstMdl.A, EstMdl.B, EstMdl.C, and EstMdl.D), and the initial state means and covariance matrix (EstMdl.Mean0 and EstMdl.Cov0), regardless of specifying Mdl explicitly. For the estimated regression coefficient, see estParams.

EstMdl does not store observed responses or predictor data. If you plan to filter (using filter), forecast (using forecast), or smooth (using smooth) using EstMdl, then you might need to supply the appropriate data.

estParams — Maximum likelihood estimates of model parametersnumeric, column vector

Maximum likelihood estimates of the model parameters known to the optimizer, returned as a numeric vector. estParams has the same dimensions as params0.

estimate arranges the estimates in estParams corresponding to unknown parameters in this order.

  1. EstMdl.A(:), that is, estimates in EstMdl.A listed column-wise

  2. EstMdl.B(:)

  3. EstMdl.C(:)

  4. EstMdl.D(:)

  5. EstMdl.Mean0

  6. EstMdl.Cov0(:)

  7. In models with predictors, estimated regression coefficients listed column-wise

Data Types: double

EstParamCov — Variance-covariance matrix of maximum likelihood estimatesmatrix

Variance-covariance matrix of maximum likelihood estimates of the model parameters known to the optimizer, returned as a matrix.

The rows and columns contain the covariances of the parameter estimates. The standard errors of the parameter estimates are the square root of the entries along the main diagonal.

estimate uses OPG to perform covariance matrix estimation.

estimate arranges the estimates in the rows and columns of EstParamCov corresponding to unknown parameters in this order.

  1. EstMdl.A(:), that is, estimates in EstMdl.A listed column-wise

  2. EstMdl.B(:)

  3. EstMdl.C(:)

  4. EstMdl.D(:)

  5. EstMdl.Mean0

  6. EstMdl.Cov0(:)

  7. In models with predictors, estimated regression coefficients listed column-wise

Data Types: double

logL — Optimized loglikelihood valuescalar

Optimized loglikelihood value, returned as a scalar.

Data Types: double

Output — Optimization informationstructure array

Optimization information, returned as a structure array.

This table describes the fields of Output.

FieldDescription
ExitFlagOptimization exit flag that describes the exit condition. For details, see fmincon and fminunc.
OptionsOptimization options that the optimizer used for numerical estimation. For details, see optimoptions.

Data Types: struct

Definitions

State-Space Model

A state-space model is a discrete-time, stochastic model that contains two sets of equations:

  • One describing how a latent process transitions in time (the state equation)

  • Another describing how an observer measures the latent process at each period (the observation equation)

Symbolically, you can write a linear, multivariate, Gaussian state-space model using the following system of equations

xt=Atxt1+BtutytZtβ=Ctxt+Dtεt,

for t = 1,...,T.

  • xt=[xt1,...,xtmt] is an mt-dimensional state vector describing the dynamics of some, possibly unobservable, phenomenon at period t.

  • yt=[yt1,...,ytnt] is an nt-dimensional observation vector describing how the states are measured by observers at period t.

  • At is the mt-by-mt – 1 state-transition matrix describing how the states at time t transition to the states at period t – 1.

  • Bt is the mt-by-kt state-disturbance-loading matrix describing how the states at period t combine with the innovations at period t.

  • Ct is the nt-by-mt measurement-sensitivity matrix describing how the observations at period t relate to the states at period t.

  • Dt is the nt-by-ht observation-innovation matrix describing how the observations at period t combine with the observation errors at period t.

  • The matrices At, Bt, Ct, and Dt are referred to as coefficient matrices, and might contain unknown parameters.

  • ut=[ut1,...,utkt] is a kt-dimensional, Gaussian, white-noise, unit-variance vector of state disturbances at period t.

  • εt=[εt1,...,εtht] is an ht-dimensional, Gaussian, white-noise, unit-variance vector of observation innovations at period t.

  • εt and ut are uncorrelated.

  • For time-invariant models,

    • Zt=[zt1zt2ztd] is row t of a T-by-d matrix of predictors Z. Each column of Z corresponds to a predictor, and each successive row to a successive period. If the observations are multivariate, then all predictors deflate each observation.

    • β is a d-by-n matrix of regression coefficients for Zt.

To write a time-invariant state-space model, drop the t subscripts of all coefficient matrices and dimensions.

Time-Invariant Model

In a time-invariant state-space model:

  • The coefficient matrices are equivalent for all periods.

  • The number of states, state disturbances, observations, and observation innovations are the same for all periods.

For example, for all t, the following system of equations

[x1,tx2,t]=[ϕ100ϕ2][x1,t1x2,t1]+[0.5002][u1,tu2,t]yt=[ϕ31][x1,tx2,t]+0.2εt

represents a time-invariant state-space model.

Time-Varying Model

In a time-varying state-space model:

  • The coefficient matrices might change from period to period.

  • The number of states, state disturbances, observations, and observation innovations might change from period to period. For example, this might happen if there is a regime shift or one of the states or observations cannot be measured during the sampling time frame. Also, you can model seasonality using time-varying models.

To illustrate a regime shift, suppose, for t = 1,..,10

[x1,tx2,t]=[ϕ100ϕ2][x1,t1x2,t1]+[0.5002][u1,tu2,t]yt=[ϕ31][x1,tx2,t]+0.2εt,

for t = 11

x1,t=[ϕ40][x1,t1x2,t1]+0.5u1,tyt=ϕ5x1,t+0.2εt,

and for t = 12,..,T

x1,t=ϕ4+0.5u1,tyt=ϕ5x1,t+0.2εt.

There are three sets of state transition matrices, whereas there are only two sets of the other coefficient matrices.

Examples

expand all

Fit a Time-Invariant State-Space Model to Data

This example generates data from a known model, and then fits a state-space model to the data.

Suppose that a latent process is this AR(1) process

$$x_t = 0.5x_{t-1} + u_t,$$

where $u_t$ is Gaussian with mean 0 and standard deviation 1.

Generate a random series of 100 observations from $x_t$, assuming that the series starts at 1.5.

T = 100;
ARMdl = arima('AR',0.5,'Constant',0,'Variance',1);
x0 = 1.5;
rng(1); % For reproducibility
x = simulate(ARMdl,T,'Y0',x0);

Suppose further that the latent process is subject to additive measurement error as indicated in the equation

$$y_t = x_t + \varepsilon_t,$$

where $\varepsilon_t$ is Gaussian with mean 0 and standard deviation 0.1.

Use the random latent state process (x) and the observation equation to generate observations.

y = x + 0.1*randn(T,1);

Together, the latent process and observation equations compose a state-space model. Supposing that the coefficients and variances are unknown parameters, the state-space model is

$$\begin{array}{c}x_t = \phi x_{t-1} + \sigma_1u_t \\ y_t = x_t +
\sigma_2\varepsilon_t.\end{array}$$

Specify the state-transition matrix. Use NaN values for unknown parameters.

A = NaN;

Specify the state-disturbance-loading coefficient matrix.

B = NaN;

Specify the measurement-sensitivity coefficient matrix.

C = 1;

Specify the observation-innovation coefficient matrix

D = NaN;

Specify the state-space model using the coefficient matrices. Also, specify the initial state mean, variance, and distribution (which is stationary).

Mean0 = 0;
Cov0 = 10;
StateType = 0;
Mdl = ssm(A,B,C,D,'Mean0',Mean0,'Cov0',Cov0,'StateType',StateType);

Mdl is an ssm model. Verify that the model is correctly specified using the display in the Command Window.

Pass the observations to estimate to estimate the parameter. Set a starting value for the parameter to params0. $\sigma_1$ and $\sigma_2$ must be positive, so set the lower bound constraints using the 'lb' name-value pair argument. Specify that the lower bound of $\phi$ is -Inf.

params0 = [0.9; 0.5; 0.1];
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf; 0; 0])
Method: Maximum likelihood (fmincon)
Sample size: 100
Logarithmic  likelihood:     -140.532
Akaike   info criterion:      287.064
Bayesian info criterion:      294.879
      |     Coeff      Std Err   t Stat    Prob  
-------------------------------------------------
 c(1) | 0.45425       0.19870   2.28611  0.02225 
 c(2) | 0.89013       0.30359   2.93205  0.00337 
 c(3) | 0.38750       0.57858   0.66975  0.50302 
      |                                          
      |  Final State   Std Dev   t Stat    Prob  
 x(1) | 1.52989       0.35621   4.29498  0.00002 

EstMdl = 


State vector length: 1
Observation vector length: 1
State disturbance vector length: 1
Observation innovation vector length: 1
Sample size supported by model: Unlimited

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equation:
x1(t) = (0.45)x1(t-1) + (0.89)u1(t)

Observation equation:
y1(t) = x1(t) + (0.39)e1(t)

Initial state distribution:

Initial state means
 x1 
  0 

Initial state covariance matrix
     x1 
 x1  10 

State types
     x1     
 Stationary 

EstMdl is an ssm model. The results of the estimation appear in the Command Window, contain the fitted state-space equations, and contain a table of parameter estimates, their standard errors, t statistics, and p-values.

You can use or display, for example the
fitted state-transition matrix using dot notation.
EstMdl.A
ans =

    0.4543

Pass EstMdl to forecast to forecast observations, or to simulate to conduct a Monte Carlo study.

Fit a Time-Varying State-Space Model to Data

This example generates data from a known model, and then fits a state-space model to the data.

Suppose that a latent process comprises an AR(2) and an MA(1) model. There are 50 periods, and the MA(1) process drops out of the model for the final 25 periods. Subsequently, the state equation for the first 25 periods is

x1,t=0.7x1,t10.2x1,t2+u1,tx2,t=u2,t+0.6u2,t1,

and for the last 25 periods, it is

x1,t=0.7x1,t10.2x1,t2+u1,t,

where u1,t and u2,t are Gaussian with mean 0 and standard deviation 1.

Generate a random series of 50 observations from x1,t and x2,t, assuming that the series starts at 1.5 and 1, respectively.

T = 50;
ARMdl = arima('AR',{0.7,-0.2},'Constant',0,'Variance',1);
MAMdl = arima('MA',0.6,'Constant',0,'Variance',1);
x0 = [1.5 1; 1.5 1];
rng(1);
x = [simulate(ARMdl,T,'Y0',x0(:,1)),...
    [simulate(MAMdl,T/2,'Y0',x0(:,2));nan(T/2,1)]];

The last 25 values for the simulated MA(1) data are NaNs.

Suppose further that the latent processes are measured using

yt=2(x1,t+x2,t)+εt,

for the first 25 periods, and

yt=2x1,t+εt

for the last 25 periods, where εt is Gaussian with mean 0 and standard deviation 1.

Use the random latent state process (x) and the observation equation to generate observations.

y = 2*nansum(x')'+randn(T,1);

Together, the latent process and observation equations compose a state-space model. Supposing that the coefficients are unknown parameters, the state-space model is

[x1,tx2,tx3,tx4,t]=[ϕ1ϕ2001000000θ10000][x1,t1x2,t1x3,t1x4,t1]+[10000101][u1,tu2,t]yt=a(x1,t+x3,t)+εt

for the first 25 periods,

[x1,tx2,t]=[ϕ1ϕ2001000][x1,t1x2,t1x3,t1x4,t1]+[10]u1,tyt=bx1,t+εt

for period 26, and

[x1,tx2,t]=[ϕ1ϕ210][x1,t1x2,t1]+[10]u1,tyt=bx1,t+εt

for the last 24 periods.

Create the function Ar2MAParamMap.m, which specifies how the parameters in params map to the state-space model matrices, the initial state values, and the type of state.

function [A,B,C,D,Mean0,Cov0,StateType] = AR2MAParamMap(params,T)
%AR2MAParamMap Time-variant state-space model parameter mapping function
%
% 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 to T/2, the
% state model is an AR(2) and an MA(1) model, and the observation model is
% the sum of the two states. From periods T/2 + 1 to T, the state model is
% just the AR(2) model.
    A1 = {[params(1) params(2) 0 0; 1 0 0 0; 0 0 0 params(3); 0 0 0 0]};
    B1 = {[1 0; 0 0; 0 1; 0 1]};
    C1 = {params(4)*[1 0 1 0]};
    Mean0 = ones(4,1);
    Cov0 = 10*eye(4);
    StateType = [0 0 0 0];
    A2 = {[params(1) params(2) 0 0; 1 0 0 0]};
    B2 = {[1; 0]};
    A3 = {[params(1) params(2); 1 0]};
    B3 = {[1; 0]};
    C3 = {params(5)*[1 0]};
    A = [repmat(A1,T/2,1);A2;repmat(A3,(T-2)/2,1)];
    B = [repmat(B1,T/2,1);B2;repmat(B3,(T-2)/2,1)];
    C = [repmat(C1,T/2,1);repmat(C3,T/2,1)];
    D = 1;
end

Specify the state-space model by passing the function AR2MAParamMap as a function handle to ssm.

Mdl = ssm(@(params)AR2MAParamMap(params,T))
Mdl = 

The state space model is implicitly defined by the function:
    @(params)AR2MAParamMap(params,T)

The software implicitly defines the state-space model. Usually, you cannot verify state-space models that you implicitly define.

Pass the observed responses (y) to estimate to estimate the parameters. Set positive, random initial values for the unknown parameters.

params0 = rand(5,1);
EstMdl = estimate(Mdl,y,params0)
Method: Maximum likelihood (fminunc)
Sample size: 50
Logarithmic  likelihood:     -114.957
Akaike   info criterion:      239.913
Bayesian info criterion:      249.473
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.47870       0.26634    1.79733  0.07229 
 c(2) |  0.00809       0.27179    0.02976  0.97626 
 c(3) |  0.55735       0.80958    0.68844  0.49118 
 c(4) |  1.62679       0.41622    3.90849  0.00009 
 c(5) |  1.90022       0.49564    3.83390  0.00013 
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | -0.81229       0.46815   -1.73511  0.08272 
 x(2) | -0.31449       0.45917   -0.68490  0.49341 

EstMdl = 


State vector length: Time-varying
Observation vector length: 1
State disturbance vector length: Time-varying
Observation innovation vector length: 1
Sample size supported by model: 50

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equations of period 1, 2, 3,..., 25:
x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t)
x2(t) = x1(t-1)
x3(t) = (0.56)x4(t-1) + u2(t)
x4(t) = u2(t)

State equations of period 26:
x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t)
x2(t) = x1(t-1)

State equations of period 27, 28, 29,..., 50:
x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t)
x2(t) = x1(t-1)


Observation equation of period 1, 2, 3,..., 25:
y1(t) = (1.63)x1(t) + (1.63)x3(t) + e1(t)

Observation equation of period 26, 27, 28,..., 50:
y1(t) = (1.90)x1(t) + e1(t)


Initial state distribution:

Initial state means
 x1  x2  x3  x4 
  1   1   1   1 

Initial state covariance matrix
     x1  x2  x3  x4 
 x1  10  0   0   0  
 x2  0   10  0   0  
 x3  0   0   10  0  
 x4  0   0   0   10 

State types
     x1          x2          x3          x4     
 Stationary  Stationary  Stationary  Stationary 

The estimated parameters are within 1 standard error of their true values, but the standard errors are quite high. Likelihood surfaces of state-space models might contain local maxima. Therefore, it is good practice to try several initial parameter values, or consider using refine.

You can pass the estimated state-space model EstMdl to forecast to forecast responses or state values.

Estimate a State-Space Model That Includes a Regression Component

Suppose that the linear relationship between the change in the unemployment rate and the nominal gross national product (nGNP) growth rate is of interest. Suppose further that the first difference of the unemployment rate is an ARMA(1,1) series. Symbolically, and in state-space form, the model is

$$\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{x_{1,t}}}\\
{{x_{2,t}}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
\phi &\theta \\
0&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_{1,t - 1}}}\\
{{x_{2,t - 1}}}
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
1&0\\
1&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{u_{1,t}}}\\
{{u_{2,t}}}
\end{array}} \right]\\
{y_t} - \beta {Z_t} = {x_{1,t}} + \sigma\varepsilon_t,
\end{array}$$

where:

  • $x_{1,t}$ is the change in the unemployment rate at time t.

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

  • $y_{1,t}$ is the observed change in the unemployment being deflated by the growth rate of nGNP ( $Z_t$).

  • $u_{1,t}$ is the Gaussian series of state disturbances having mean 0 and standard deviation 1.

  • $\varepsilon_t$ is the Gaussian series of observation innovations having mean 0 and standard deviation $\sigma$.

Load the Nelson-Plosser data set, which contains the unemployment rate and nGNP series, among other things.

load Data_NelsonPlosser

Preprocess the data by taking the natural logarithm of the nGNP series, and the first difference of each. Also, remove the starting NaN values from each series.

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
u = DataTable.UR(~isNaN);
T = size(gnpn,1);                          % Sample size
Z = [ones(T-1,1) diff(log(gnpn))];
y = diff(u);

This example proceeds using series without NaN values. However, using the Kalman filter framework, the software can accommodate series containing missing values.

Specify the state-transition coefficient matrix.

A = [NaN NaN; 0 0];

Specify the state-disturbance-loading coefficient matrix.

B = [1 0;1 0];

Specify the measurement-sensitivity coefficient matrix.

C = [1 0];

Specify the observation-innovation coefficient matrix.

D = NaN;

Specify the state-space model using ssm.

Mdl = ssm(A,B,C,D);

Estimate the model parameters. Specify the regression component and its initial value for optimization using the 'Predictors' and 'Beta0' name-value pair arguments, respectively. Display the estimates and all optimization diagnostic information. Restrict the estimate of $\sigma$ to all positive, real numbers.

params0 = [0.3 0.2 0.1]; % Chosen arbitrarily
EstMdl = estimate(Mdl,y,params0,'Predictors',Z,'Display','full',...
    'Beta0',randn(2,1),'lb',[-Inf,-Inf,0,-Inf,-Inf]);
____________________________________________________________
   Diagnostic Information

Number of variables: 5

Functions 
Objective:                            @(c)-fML(c,Mdl,Y,Predictors,unitFlag,sqrtFlag,mexFlag,mexTvFlag,tol,ind)
Gradient:                             finite-differencing
Hessian:                              finite-differencing (or Quasi-Newton)

Constraints
Nonlinear constraints:                do not exist
 
Number of linear inequality constraints:    0
Number of linear equality constraints:      0
Number of lower bound constraints:          1
Number of upper bound constraints:          0

Algorithm selected
   interior-point


____________________________________________________________
   End diagnostic information
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       6    2.815114e+02    0.000e+00    4.850e+01
    1      20    2.792263e+02    0.000e+00    4.119e+01    1.387e-01
    2      27    2.722399e+02    0.000e+00    4.692e+01    2.302e-01
    3      34    2.425661e+02    0.000e+00    1.784e+01    1.268e+01
    4      40    1.919274e+02    0.000e+00    2.966e+01    5.524e+00
    5      47    1.505801e+02    0.000e+00    1.120e+02    6.986e+00
    6      56    1.344070e+02    0.000e+00    2.316e+02    8.287e-01
    7      69    1.305133e+02    0.000e+00    1.592e+02    2.095e-01
    8      76    1.281543e+02    0.000e+00    1.018e+02    2.068e-01
    9      83    1.220398e+02    0.000e+00    7.499e+01    3.084e+00
   10      89    1.174846e+02    0.000e+00    4.571e+01    2.133e+00
   11      95    1.092141e+02    0.000e+00    1.012e+01    2.238e+00
   12     101    1.089658e+02    0.000e+00    4.504e+00    7.843e-01
   13     107    1.085915e+02    0.000e+00    2.260e+01    5.076e-01
   14     119    1.075450e+02    0.000e+00    3.910e+01    2.497e-01
   15     144    1.075445e+02    0.000e+00    1.397e+01    4.118e-05
   16     151    1.075443e+02    0.000e+00    1.402e+01    8.240e-05
   17     160    1.073947e+02    0.000e+00    1.997e+02    3.142e-01
   18     168    1.072711e+02    0.000e+00    7.103e+02    1.675e-01
   19     180    1.072652e+02    0.000e+00    2.344e+02    2.985e-05
   20     187    1.072652e+02    0.000e+00    2.345e+02    5.931e-05
   21     195    1.071561e+02    0.000e+00    1.586e+03    7.935e-02
   22     202    1.070592e+02    0.000e+00    4.680e+03    1.377e-01
   23     213    1.070398e+02    0.000e+00    3.341e+02    7.989e-06
   24     220    1.070398e+02    0.000e+00    3.354e+02    1.507e-05
   25     226    1.069775e+02    0.000e+00    5.148e+04    2.493e-01
   26     238    1.069581e+02    0.000e+00    2.070e+04    2.576e-06
   27     245    1.069581e+02    0.000e+00    2.070e+04    2.373e-06
   28     251    1.068345e+02    0.000e+00    1.393e+04    1.237e-01
   29     259    1.064422e+02    0.000e+00    3.871e+04    3.056e-01
   30     269    1.064370e+02    0.000e+00    1.421e+04    5.824e-07

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31     277    1.061184e+02    0.000e+00    1.088e+06    3.056e-01
   32     293    1.057993e+02    0.000e+00    2.644e+05    1.870e-07

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.



Warning: Covariance matrix of estimators cannot be computed precisely due to
inversion difficulty. Check parameter identifiability. Also try different
starting values and other options to compute the covariance matrix. 

Method: Maximum likelihood (fmincon)
Sample size: 61
Logarithmic  likelihood:     -105.799
Akaike   info criterion:      221.599
Bayesian info criterion:      232.153
           |      Coeff       Std Err       t Stat         Prob  
-----------------------------------------------------------------
 c(1)      |   1.00000       0.00000        4.32645e+06   0      
 c(2)      |  -1.00007       0.00005   -20140.44821       0      
 c(3)      |   0.84579       0.10647        7.94389       0      
 y <- z(1) |   1.03050       0.24444        4.21584      0.00002 
 y <- z(2) | -24.54687       1.42559      -17.21879       0      
           |                                                     
           |    Final State   Std Dev        t Stat        Prob  
 x(1)      |   0.98612       0.64842        1.52081      0.12831 
 x(2)      |   0.76116       0.65092        1.16935      0.24227 

Optimization information and a table of estimates and statistics output to the Command Window. EstMdl is an ssm model, and you can access its properties using dot notation.

Compare Estimates from State-Space Model Filtering Methods

The software implements the Kalman filter using the covariance filter by default, but you can specify to use the square-root filter instead. This example compares estimates from each method using simulated data.

Suppose that a latent process is an AR(1). Subsequently, the state equation is

$$x_t = 0.5x_{t-1} + u_t,$$

where $u_t$ is Gaussian with mean 0 and standard deviation 0.3.

Generate a random series of 100 observations from $x_t$, assuming that the series starts at 1.5.

T = 100;
ARMdl = arima('AR',0.5,'Constant',0,'Variance',0.3^2);
x0 = 1.5;
rng(1); % For reproducibility
x = simulate(ARMdl,T,'Y0',x0);

Suppose further that the latent process is subject to additive measurement error. Subsequently, the observation equation is

$$y_t = x_t + \varepsilon_t,$$

where $\varepsilon_t$ is Gaussian with mean 0 and standard deviation 0.1.

Use the random latent state process (x) and the observation equation to generate observations.

y = x + 0.1*randn(T,1);

Together, the latent process and observation equations compose a state-space model. Supposing that the coefficients and variances are unknown parameter, the state-space model is

$$\begin{array}{c}x_t = \phi x_{t-1} + \sigma_1u_t \\ y_t = x_t +&#xA;\sigma_2\varepsilon_t\end{array}$$

Specify the state-transition coefficient matrix. Use NaN values for unknown parameters.

A = NaN;

Specify the state-disturbance-loading coefficient matrix.

B = NaN;

Specify the measurement-sensitivity coefficient matrix.

C = 1;

Specify the observation-innovation coefficient matrix.

D = NaN;

Specify the state-space model using the coefficient matrices. Also, specify the initial state mean, variance, and distribution (which is stationary).

Mean0 = 0;
Cov0 = 10;
StateType = 0;
Mdl = ssm(A,B,C,D,'Mean0',Mean0,'Cov0',Cov0,'StateType',StateType);

Mdl is an ssm model.

Estimate the parameters using estimate two ways:

  • Using the default, simple Kalman filter

  • Using the square root filter variation

In both cases, specify that no output should be returned to the Command Window. This is good practice if you plan on running estimate multiple times (such as a Monte Carlo simulation).

params0 = [10,10,10];
[~,estParamsSKF,EstParamCovSKF,logLSKF,OutputSKF] = estimate(Mdl,y,params0,...
    'Display','off');
[~,estParamsSR,EstParamCovSR,logLSR,OutputSR] = estimate(Mdl,y,params0,...
    'Squareroot',true,'Display','off');

Check that the algorithms converged properly by printing the exit flag properties of OutputSKF and OutputSR.

exitFlagSKF = OutputSKF.ExitFlag
exitFlagSR = OutputSR.ExitFlag
exitFlagSKF =

     1


exitFlagSR =

     1

Both algorithms have an exit flag of 1, which indicates that the software met the convergence criteria.

Compare the estimates from each algorithm.

fprintf('\n Parameter Estimates\n')
table(estParamsSKF',estParamsSR','VariableNames',...
    {'SimpleKalmanFilter','SquarerootFilter'})
fprintf('\nEstimated Parameter Covariance Matrix\n')
table(EstParamCovSKF,EstParamCovSR,'VariableNames',...
    {'SimpleKalmanFilter','SquarerootFilter'})
 Parameter Estimates

ans = 

    SimpleKalmanFilter    SquarerootFilter
    __________________    ________________

     0.51057               0.51057        
     0.23436               0.23436        
    -0.17904              -0.17904        


Estimated Parameter Covariance Matrix

ans = 

            SimpleKalmanFilter                      SquarerootFilter          
    ___________________________________    ___________________________________

     0.036669    -0.013302    -0.014012     0.036669    -0.013302    -0.014012
    -0.013302    0.0070187    0.0072533    -0.013302    0.0070187    0.0072533
    -0.014012    0.0072533    0.0089019    -0.014012    0.0072533    0.0089019

In this case, the results are the same.

If you use the default, covariance filter method, and you run into numerical problems during estimation, filtering, or smoothing, try using the squareroot method.

Algorithms

  • For explicitly defined state-space models, the software applies all predictors to each response series. However, each response series has its own set of regression coefficients.

  • The software estimates regression coefficients along with all other state-space model parameters. The software is flexible enough to allow applying constraints to the regression coefficients using the constrained optimization options for fminunc.

  • The software passes the name-value pair arguments Options, Aineq, bineq, Aeq, beq, lb, and ub directly to the optimizer fmincon or fminunc.

  • If you do not specify optimization constraints, then the software uses fminunc for unconstrained numerical estimation. If you specify any pair of optimization constraints, then the software uses fmincon for constrained numerical estimation. For either type of optimization, optimization options you set using the name-value pair argument Options must be consistent with the options of the optimization algorithm.

  • If you set 'Univariate',true, then, during the filtering algorithm, the software sequentially updates rather then updating all at once. This might accelerate parameter estimation, especially for a low-dimensional, time-invariant model.

  • Suppose that you want to create a state-space model using a parameter-to-matrix mapping function with this signature

    [A,B,C,D,Mean0,Cov0,StateType,DeflateY] = paramMap(params,Y,Z)

    and you specify the model using an anonymous function

    Mdl = ssm(@(params)paramMap(params,Y,Z))

    The observed responses Y and predictor data Z are not input arguments in the anonymous function. If Y and Z exist in the MATLAB Workspace before creating Mdl, then the software establishes a link to them. Otherwise, if you pass Mdl to estimate, the software throws an error.

    The link to the data established by the anonymous function overrides all other corresponding input argument values of estimate. This distinction is important particularly when conducting a rolling window analysis. For details, see Rolling-Window Analysis of Time-Series Models.

References

[1] Durbin J., and S. J. Koopman. Time Series Analysis by State Space Methods. 2nd ed. Oxford: Oxford University Press, 2012.

Was this topic helpful?