filter

Suppose that a latent process is an AR(1). The state equation is

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

where $${\epsilon }_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)

Mdl = 
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.50)x1(t-1) + u1(t)

Observation equation:
y1(t) = x1(t) + (0.75)e1(t)

Initial state distribution:

Initial state means
 x1 
  0 

Initial state covariance matrix
     x1   
 x1  1.33 

State types
     x1     
 Stationary 

Mdl is an ssm model. Verify that the model is correctly specified using the display in the Command Window. The software infers that the state process is stationary. Subsequently, the software sets the initial state mean and covariance to the mean and variance of the stationary distribution of an AR(1) model.

Filter states for periods 1 through 100. Plot the true state values and the filtered state estimates.

filteredX = filter(Mdl,y);

figure
plot(1:T,x,'-k',1:T,filteredX,':r','LineWidth',2)
title({'State Values'})
xlabel('Period')
ylabel('State')
legend({'True state values','Filtered state values'})

The true values and filter estimates are approximately the same.

Suppose that the linear relationship between the change in the unemployment rate and the nominal gross national product (nGNP) growth rate is of interest. Suppose further that the first difference of the unemployment rate is an ARMA(1,1) series. Symbolically, and in state-space form, the model is

where:

Load the Nelson-Plosser data set, which contains the unemployment rate and nGNP series, among other things.

load Data_NelsonPlosser

Preprocess the data by taking the natural logarithm of the nGNP series, and the first difference of each series. Also, remove the starting NaN values from each series.

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
u = DataTable.UR(~isNaN);
T = size(gnpn,1);                          % Sample size
Z = [ones(T-1,1) diff(log(gnpn))];
y = diff(u);

Though this example removes missing values, the software can accommodate series containing missing values in the Kalman filter framework.

Specify the coefficient matrices.

A = [NaN NaN; 0 0];
B = [1; 1];
C = [1 0];
D = NaN;

Specify the state-space model using ssm.

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

Estimate the model parameters, and use a random set of initial parameter values for optimization. Specify the regression component and its initial value for optimization using the 'Predictors' and 'Beta0' name-value pair arguments, respectively. Restrict the estimate of $$\sigma$$ to all positive, real numbers.

params0 = [0.3 0.2 0.2];
[EstMdl,estParams] = estimate(Mdl,y,params0,'Predictors',Z,...
    'Beta0',[0.1 0.2],'lb',[-Inf,-Inf,0,-Inf,-Inf]);

Method: Maximum likelihood (fmincon)
Sample size: 61
Logarithmic  likelihood:     -99.7245
Akaike   info criterion:      209.449
Bayesian info criterion:      220.003
           |      Coeff       Std Err    t Stat     Prob  
----------------------------------------------------------
 c(1)      |  -0.34098       0.29608    -1.15164  0.24948 
 c(2)      |   1.05003       0.41377     2.53771  0.01116 
 c(3)      |   0.48592       0.36790     1.32079  0.18657 
 y <- z(1) |   1.36121       0.22338     6.09358   0      
 y <- z(2) | -24.46711       1.60018   -15.29024   0      
           |                                              
           |    Final State   Std Dev     t Stat    Prob  
 x(1)      |   1.01264       0.44690     2.26592  0.02346 
 x(2)      |   0.77718       0.58917     1.31912  0.18713 

EstMdl is an ssm model, and you can access its properties using dot notation.

Filter the estimated state-space model. EstMdl does not store the data or the regression coefficients, so you must pass in them in using the name-value pair arguments 'Predictors' and 'Beta', respectively. Plot the estimated, filtered states. Recall that the first state is the change in the unemployment rate, and the second state helps build the first.

filteredX = filter(EstMdl,y,'Predictors',Z,'Beta',estParams(end-1:end));

figure
plot(dates(end-(T-1)+1:end),filteredX(:,1));
xlabel('Period')
ylabel('Change in the unemployment rate')
title('Filtered Change in the Unemployment Rate')

Consider nowcasting the model in Filter States of Time-Invariant State-Space Model.

Generate a random series of 100 observations from $$x_t$$.

T = 100;
A = 0.5;
B = 1;
C = 1;
D = 0.75;
Mdl = ssm(A,B,C,D);
rng(1); % For reproducibility
y = simulate(Mdl,T);

Suppose the final 10 observations are in the forecast horizon.

fh = 10;
yf = y((end-fh+1):end); % Holdout sample responses
y = y(1:end-fh);        % In-sample responses

Filter the observations through the model to obtain filtered states for each period.

[xhat,logL,output] = filter(Mdl,y);
xhatvar = output.FilteredStatesCov;

xhat and xhatvar are 90-by-1 vectors of in-sample filtered states and corresponding variances, respectively. xhat(t) is the estimate of $\mathit{E}\left({\mathit{x}}_{\mathit{t}}|{\mathit{y}}_1,...,{\mathit{y}}_{\mathit{t}}\right)$, and xhatvar is the estimate of $\mathrm{Var}\left({\mathit{x}}_{\mathit{t}}|{\mathit{y}}_1,...,{\mathit{y}}_{\mathit{t}}\right)$.

Call filter again, but specify the real-time update option.

[xhatRT,logLRT,outputRT] = filter(Mdl,y,'RealTimeUpdate',true);
xhatRTvar = outputRT.FilteredStatesCov;

xhatRT and xhatRTvar are scalars representing the estimate of $\mathit{E}\left({\mathit{x}}_{90}|{\mathit{y}}_1,...,{\mathit{y}}_{90}\right)$ and its corresponding variance, respectively.

Compare the filtered states and variances of period 90.

tol = 1e-10;
areMeansEqual = (xhat(end) - xhatRT) < tol

areMeansEqual = logical
   1

areVarsEqual = (xhatvar(end) - xhatRTvar) < tol

areVarsEqual = logical
   0

areLogLsEqual = (logL - logLRT) < tol;

In the last period, the filtered states, their variances, and the loglikelihoods are equal.

Nowcast the model into the forecast horizon by performing this procedure for each successive period:

% Initialize state and variance
state0 = xhatRT;        
var0 = xhatRTvar;

% Preallocate
xhatRTF = zeros(fh,1);
xhatRTvarF = zeros(fh,1);

for j = 1:fh
    Mdl.Mean0 = state0;
    Mdl.Cov0 = var0;
    [xhatRTF(j),~,outputRT] = filter(Mdl,yf(j),'RealTimeUpdate',true); % Alternatively, use update
    xhatRTvarF(j) = outputRT.FilteredStatesCov;
    state0 = xhatRTF(j);
    var0 = xhatRTvarF(j);
end

Plot the data and nowcasts.

figure
plot((T-fh-20):T,[y(end-20:end); yf],'b-',(T-fh+1):T,xhatRTF,'r*-')
legend(["Data" "Nowcasts"],'Location',"best")