%% FEVD of Estimated Model
% Simulate data from a known model, fit the data to a state-space model, and then estimate the FEVD of the measurement variables.
% Consider the time series decomposition y_t =\tau_t +c_t, where \tau_t is a random walk with drift representing the trend component, and c_t is an AR(1) model representing the cycle component.
% \tau_t=3 + \tau_{t-1}+u_{1,t}\\
% c_t=0.5c_{t-1}+2u_{2,t}.
% The model in state-space notation is
% \begin{array}{l}
% {x_t} = \left[ {\begin{array}{*{20}{c}}
% {{\tau _t}}\\
% {{d_t}}\\
% {{c_t}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% 1&3&0\\
% 0&1&0\\
% 0&0&{0.5}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{\tau _{t - 1}}}\\
% {{d_{t - 1}}}\\
% {{c_{t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1&0\\
% 0&0\\
% 0&2
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{u_{1,t}}}\\
% {{u_{2,t}}}
% \end{array}} \right]\\
% {y_t} = \left[ {\begin{array}{*{20}{c}}
% 1&0&1
% \end{array}} \right]{x_t},\end{array}
% where d_t is a dummy state representing the drift parameter, which is 1 for all t.
% Simulate 500 observations from the true model.

rng(1); % For reproducibility
ADGP = [1 3 0; 0 1 0; 0 0 0.5];
BDGP = [1 0; 0 0; 0 2];
CDGP = [1 0 1];
DGP = ssm(ADGP,BDGP,CDGP,'StateType',[2 1 0]);
y = simulate(DGP,500);

% Assume that the drift constant, disturbance variances, and AR coefficient are unknown. Explicitly create a state-space model template for estimation that represents the model by replacing the unknown parameters in the model with NaN.

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

% Fit the model template to the data. Specify a set of positive, random standard Gaussian starting values for the four model parameters. Return the estimated model and vector of parameter estimates.

[EstMdl,estParams] = estimate(Mdl,y,abs(randn(4,1)),'Display','off')

% EstMdl is a fully specified ssm model object. Model estimates are close to their true values.
% Compute and plot the FEVD of the measurement variable. Specify the model template Mdl and the vector of estimated parameters estParams.

Decomposition = fevd(Mdl,'Params',estParams);

bar(Decomposition,'stacked')
xlabel('Period')
ylabel('Variance decompositions of $y_{t}$','Interpreter','latex')
legend('$u_{1,t}$','$u_{2,t}$','Interpreter','latex')

% Noise in the cyclical component dominates the volatility of the measurement variable in low lags, with increasing contribution from the trend component noise as the lag increases.
% Copyright 2021 The MathWorks, Inc.