%% Estimate Posterior Distribution of States in State-Space Model
% The simsmooth function draws random samples from the distribution of smoothed states, or the distribution of a state given all of the data and parameters. This is the definition of posterior distribution of a state. Suppose that a latent process is an AR(1). The state equation is
% 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. The observation equation is
% y_t = x_t + \varepsilon_t,
% where \varepsilon_t is Gaussian with mean 0 and standard deviation 0.75. Together, the latent process and observation equations compose a state-space model.
%
% Use the random latent state process (x) and the observation equation to generate observations.

y = x + 0.75*randn(T,1);

%
% Specify the four coefficient matrices.

A = 0.5;
B = 1;
C = 1;
D = 0.75;

%
% Specify the state-space model using the coefficient matrices.

Mdl = ssm(A,B,C,D);

%
% Smooth the states of the state space model.

xsmooth = smooth(Mdl,y);

%
% Draw 1000 paths from the posterior distribution of x_1.

N = 1000;
SimX = simsmooth(Mdl,y,'NumPaths',N);

% SimX is a 100-by- 1-by- 1000 array. Rows correspond to periods, columns correspond to individual states, and leaves correspond to separate paths.
%
% Because SimX has a singleton dimension, collapse it so that its leaves correspond to the columns using squeeze.

SimX = squeeze(SimX);

%
% Compute the mean, standard deviation, and 95% confidence intervals of the state at each period.

xbar = mean(SimX,2);
xstd = std(SimX,[],2);
ci = [xbar - 1.96*xstd, xbar + 1.96*xstd];

%
% Plot the smoothed states, and the means and 95% confidence intervals of the draws at each period.

figure;
plot(xsmooth,'k','LineWidth',2);
hold on;
plot(xbar,'--r','LineWidth',2);
plot(1:T,ci(:,1),'--r',1:T,ci(:,2),'--r');
legend('Smoothed states','Simulation Mean','95% CIs');
title('Smooth States and Simulation Statistics');
xlabel('Period')

% Copyright 2012 The MathWorks, Inc.