%% Temporal Covariances of Estimated Model
% Simulate data from a known model, fit a model to the data, and then compare sample and model-implied covariances.
%
%% Simulate Data
% Explicitly create the state-space model
% \begin{array}{*{20}{l}}
% {{x_{1,t}}} = 0.9x_{1,t-1} + 0.2u_{1,t}\\
% {{x_{2,t}}} = 0.1x_{1,t-1}-0.3x_{2,t-1} + u_{2,t}\\
% {{y_{1,t}}} = x_{1,t}  +\varepsilon_{1,t}\\
% {{y_{2,t}}} = x_{1,t} + x_{2,t} +\varepsilon_{2,t}.\\
% \end{array}}

ADGP = [0.9 0; 0.1 -0.3];
BDGP = [0.2 0; 0 1];
CDGP = [1 0; 1 1];
DDGP = eye(2);
DGP = ssm(ADGP,BDGP,CDGP,DDGP,'StateType',[0 0]);

% Generate a random path of measurements of length 500 from the model.

rng(1); % For reproducibility
numobs = 500;
Y = simulate(DGP,numobs);

%
%% Fit Model to Data
% Create a state-space model template to fit to the data by replacing each nonzero state parameter of the data-generating process with a NaN value.

A = [NaN 0; NaN NaN];
B = [NaN 0; 0 NaN];
Mdl = ssm(A,B,CDGP,DDGP,'StateType',[0  0]);

% Fit the model template to the data. Specify a random set of positive starting values. Return the vector of estimated parameters.

[~,estParams] = estimate(Mdl,Y,abs(rand(5,1)));

%
%% Compute Covariances
% Compute model-implied temporal covariances of the measurement variables by passing the state-space model template and estimated parameters to corr. Return the covariances instead of the correlations.

Covyy = corr(Mdl,'Params',estParams,'Covariance',true)

% Covyy is a 2-by-2-by-2 array of temporal covariances of the measurement variables.
% Compute the sample covariances of the measurement variables and their first lags.

AugData = lagmatrix(Y,[0 1]);
SCovyy = cov(AugData(2:end,:));

%
%% Compare Covariances
% Compare the model-implied temporal covariances and the sample covariances.

names = ["y_1t" "y_2t" "y_1t-1" "y_2t-1"];

Covy1y = squeeze(Covyy(:,1,:))';
Covy2y = squeeze(Covyy(:,2,:))';
CovyLag1 = [Covy1y(:,2) Covy2y(:,2) Covy1y(:,1) Covy2y(:,1)]';
ModelCovariances = array2table([Covy1y(:) Covy2y(:) CovyLag1],'RowNames',names,...
    'VariableNames',names)
SampleCovariances = array2table(SCovyy,'RowNames',names,'VariableNames',names)

% The model-implied and sample covariances appear to be similar in magnitude. Note that the covariances are invariant to the reference time; for example, Cov(y_{1,t} ,y_{2,t}) = Cov(y_{1,t-1} ,y_{2,t-1}).
% Copyright 2021 The MathWorks, Inc.