Estimate Time-Invariant State-Space Model

This example shows how to generate data from a known model, specify a state-space model containing unknown parameters corresponding to the data generating process, and then fit the state-space model to the data.

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

${x}_{t}=0.5{x}_{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}+{\epsilon }_{t},$

where ${\epsilon }_{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}=\varphi {x}_{t-1}+{\sigma }_{1}{u}_{t}\\ {y}_{t}={x}_{t}+{\sigma }_{2}{\epsilon }_{t}.\end{array}$

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

A = NaN;

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 $\varphi$ 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-space model type: ssm

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.