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.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+{\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;

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 $$\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.52990       0.35620   4.29499  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.