%% Smooth Time-Varying State-Space Model
% This example shows how to generate data from a known model, fit a
% state-space model to the data, and then smooth the states.
%%
% 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
%
% $$\begin{array}{l}
% {x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}}\\
% {x_{2,t}} = {u_{2,t}} + 0.6{u_{2,t - 1}},
% \end{array}$$
%
% and for the last 25 periods, it is 
%
% $${x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}},$$
%
% where $u_{1,t}$ and $u_{2,t}$ are Gaussian with mean 0 and standard
% deviation 1.
%%
% Assuming that the series starts at 1.5 and 1, respectively, generate a
% random series of 50 observations from $x_{1,t}$ and $x_{2,t}$.

% Copyright 2015 The MathWorks, Inc.

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 |NaN| values.
%%
% Suppose further that the latent processes are measured using
%
% $${y_t} = 2\left( {{x_{1,t}} + {x_{2,t}}} \right) + {\varepsilon _t},$$
%
% for the first 25 periods, and
%
% $${y_t} = 2{x_{1,t}} + {\varepsilon _t}$$
%
% for the last 25 periods, where $\varepsilon_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
% 
% $$\begin{array}{l}
% \left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}\\
% {{x_{3,t}}}\\
% {{x_{4,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% {{\phi _1}}&{{\phi _2}}&0&0\\
% 1&0&0&0\\
% 0&0&0&{{\theta _1}}\\
% 0&0&0&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t - 1}}}\\
% {{x_{2,t - 1}}}\\
% {{x_{3,t - 1}}}\\
% {{x_{4,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1&0\\
% 0&0\\
% 0&1\\
% 0&1
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{u_{1,t}}}\\
% {{u_{2,t}}}
% \end{array}} \right]\\
% {y_t} = a({x_{1,t}} + {x_{3,t}}) + {\varepsilon _t}
% \end{array}$$
%
% for the first 25 periods,
%
% $$\begin{array}{l}
% \left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% {{\phi _1}}&{{\phi _2}}&0&0\\
% 1&0&0&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t - 1}}}\\
% {{x_{2,t - 1}}}\\
% {{x_{3,t - 1}}}\\
% {{x_{4,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1\\
% 0
% \end{array}} \right]{u_{1,t}}\\
% {y_t} = b{x_{1,t}} + {\varepsilon _t}
% \end{array}$$
%
% for period 26, and
%
% $$\begin{array}{l}
% \left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% {{\phi _1}}&{{\phi _2}}\\
% 1&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
% \end{array}} \right]{u_{1,t}}\\
% {y_t} = b{x_{1,t}} + {\varepsilon _t}
% \end{array}$$
%
% for the last 24 periods.
%%
% Write a function that specifies how the parameters in |params| map to the
% state-space model matrices, the initial state values, and the type of
% state.
%
% <include>AR2MAParamMap.m</include>
%
%%
% Save this code as a file named |AR2MAParamMap| on your MATLAB(R) path.
%%
% Create the state-space model by passing the function |AR2MAParamMap| as a
% function handle to |ssm|.
Mdl = ssm(@(params)AR2MAParamMap(params,T));
%%
% |ssm| implicitly creates the state-space model. Usually, you
% cannot verify an implicitly defined state-space model.
%%
% Pass the observed responses (|y|) to |estimate| to estimate the
% parameters. Specify an arbitrary set of positive initial values for the
% unknown parameters.
params0 = 0.1*ones(5,1);
EstMdl = estimate(Mdl,y,params0);
%%
% |EstMdl| is an |ssm| model containing the estimated coefficients.
% 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|.
%%
% Smooth the states and estimate the variance-covariance matrices of the
% smoothed states by passing |EstMdl| and the observed responses to
% |smooth|.
[~,~,Output]= smooth(EstMdl,y);
%%
% |Output| is a |T|-by-1 structure array containing the smoothed states and
% their variance-covariance matrices, among other things.
%%
% Extract the smoothed states and their variances from the cell arrays.
% Recall that the two, different states are in positions 1 and 3. The
% states in positions 2 and 4 help specify the processes of interest.
stateIndx = [1 3]; % State Indices of interest
SmoothedStates = NaN(T,numel(stateIndx));
SmoothedStatesCov = NaN(T,numel(stateIndx));

for t = 1:T
    maxInd1 = size(Output(t).SmoothedStates,1);
		maxInd2 = size(Output(t).SmoothedStatesCov,1);
    mask1 = stateIndx <= maxInd1;
		mask2 = stateIndx <= maxInd2;
    SmoothedStates(t,mask1) = ...
        Output(t).SmoothedStates(stateIndx(mask1),1);
		SmoothedStatesCov(t,mask2) = ...
        diag(Output(t).SmoothedStatesCov(stateIndx(mask2),...
        stateIndx(mask2)));
end
%%
% Plot the true state values, the smoothed state values, and their
% individual 95% Wald-type confidence intervals for each model.
AR2SSCIlb = SmoothedStates(:,1) - 1.95*sqrt(SmoothedStatesCov(:,1));
AR2SSCIub = SmoothedStates(:,1) + 1.95*sqrt(SmoothedStatesCov(:,1));
AR2SSIntervals = [AR2SSCIlb AR2SSCIub];

MA1SSCIlb = SmoothedStates(:,2) - 1.95*sqrt(SmoothedStatesCov(:,2));
MA1SSCIub = SmoothedStates(:,2) + 1.95*sqrt(SmoothedStatesCov(:,2));
MA1SSIntervals = [MA1SSCIlb MA1SSCIub];

figure
plot(1:T,x(:,1),'-k',1:T,SmoothedStates(:,1),':r',...
		1:T,AR2SSIntervals,'--b','LineWidth',2);
title('AR(2) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Smoothed state values',...
    '95% Confidence Intervals'});

figure
plot(1:T,x(:,2),'-k',1:T,SmoothedStates(:,2),':r',...
		1:T,MA1SSIntervals,'--b','LineWidth',2);
title('MA(1) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Smoothed state values',...
    '95% Confidence Intervals'});
