smooth

Suppose that a latent process is a random walk. 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;
x0 = 1.5;
rng(1); % For reproducibility
u = randn(T,1);
x = cumsum([x0;u]);
x = x(2:end);

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 = 1;
B = 1;
C = 1;
D = 0.75;

Create the diffuse state-space model using the coefficient matrices. Specify that the initial state distribution is diffuse.

Mdl = dssm(A,B,C,D,'StateType',2)

Mdl = 
State-space model type: dssm

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) = 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  Inf 

State types
    x1   
 Diffuse 

Mdl is an dssm model. Verify that the model is correctly specified using the display in the Command Window.

Smooth states for periods 1 through 100. Plot the true state values and the smoothed state estimates.

SmoothedX = smooth(Mdl,y);

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

The true values and smoothed estimates are approximately the same.

Suppose that the linear relationship between unemployment rate and the nominal gross national product (nGNP) is of interest. Suppose further that unemployment rate is an AR(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 removing the starting NaN values from each series.

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
y = diff(DataTable.UR(~isNaN));
T = size(gnpn,1);                          % The sample size
Z = price2ret(gnpn);

This example continues using the series without NaN values. However, using the Kalman filter framework, the software can accommodate series containing missing values.

Specify the coefficient matrices.

A = NaN;
B = NaN;
C = 1;

Create the state-space model using dssm by supplying the coefficient matrices and specifying that the state values come from a diffuse distribution. The diffuse specification indicates complete ignorance about the moments of the initial distribution.

StateType = 2;
Mdl = dssm(A,B,C,'StateType',StateType);

Estimate the parameters. Specify the regression component and its initial value for optimization using the 'Predictors' and 'Beta0' name-value pair arguments, respectively. Display the estimates and all optimization diagnostic information. Restrict the estimate of $$\sigma$$ to all positive, real numbers.

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

Method: Maximum likelihood (fmincon)
Effective Sample size:             60
Logarithmic  likelihood:     -110.477
Akaike   info criterion:      226.954
Bayesian info criterion:      233.287
           |      Coeff       Std Err    t Stat    Prob 
--------------------------------------------------------
 c(1)      |   0.59436       0.09408     6.31738  0     
 c(2)      |   1.52554       0.10758    14.17991  0     
 y <- z(1) | -24.26161       1.55730   -15.57930  0     
           |                                            
           |    Final State   Std Dev     t Stat   Prob 
 x(1)      |   2.54764        0           Inf     0     

EstMdl is a dssm model, and you can access its properties using dot notation.

Smooth the estimated diffuse 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 smoothed states.

SmoothedX = smooth(EstMdl,y,'Predictors',Z,'Beta',estParams(end));

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

Estimate a diffuse state-space model, smooth the states, and then extract other estimates from the Output output argument.

Consider the diffuse state-space model

The state variable $$x_{1,t}$$ is an AR(1) model with autoregressive coefficient $$\varphi$$. $$x_{2,t}$$ is a random walk. The disturbances $$u_{1,t}$$ and $$u_{2,t}$$ are independent Gaussian random variables with mean 0 and standard deviations $${\sigma }_1$$ and $${\sigma }_2$$, respectively. The observation $$y_t$$ is the error-free sum of $$x_{1,t}$$ and $$x_{2,t}$$.

Generate data from the state-space model. To simulate the data, suppose that the sample size $$T=100$$, $$\varphi =0.6$$, $${\sigma }_1=0.2$$, $${\sigma }_2=0.1$$, and $$x_{1,0}=x_{2,0}=2$$.

rng(1); % For reproducibility
T = 100;
ARMdl = arima('AR',0.6,'Constant',0,'Variance',0.2^2);
x1 = simulate(ARMdl,T,'Y0',2);
u3 = 0.1*randn(T,1);
x3 = cumsum([2;u3]);
x3 = x3(2:end);
y = x1 + x3;

Specify the coefficient matrices of the state-space model. To indicate unknown parameters, use NaN values.

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

Create a diffuse state-space model that describes the model above. Specify that $$x_{1,t}$$ and $$x_{2,t}$$ have diffuse initial state distributions.

StateType = [2 2];
Mdl = dssm(A,B,C,'StateType',StateType);

Estimate the unknown parameters of Mdl. Choose initial parameter values for optimization. Specify that the standard deviations are constrained to be positive, but all other parameters are unconstrained using the 'lb' name-value pair argument.

params0 = [0.01 0.1 0.01]; % Initial values chosen arbitrarily
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf 0 0]);

Method: Maximum likelihood (fmincon)
Effective Sample size:             98
Logarithmic  likelihood:      3.44283
Akaike   info criterion:    -0.885655
Bayesian info criterion:      6.92986
      |     Coeff      Std Err   t Stat     Prob  
--------------------------------------------------
 c(1) | 0.54134       0.20494    2.64145  0.00826 
 c(2) | 0.18439       0.03305    5.57897   0      
 c(3) | 0.11783       0.04347    2.71039  0.00672 
      |                                           
      |  Final State   Std Dev    t Stat    Prob  
 x(1) | 0.24884       0.17168    1.44943  0.14722 
 x(2) | 1.73762       0.17168   10.12121   0      

The parameters are close to their true values.

Smooth the states of EstMdl, and request all other available output.

[X,logL,Output] = smooth(EstMdl,y);

X is a T-by-2 matrix of smoothed states, logL is the final optimized log-likelihood value, and Output is a structure array containing various estimates that the Kalman filter requires. List the fields of output using fields.

fields(Output)

ans = 9x1 cell
    {'LogLikelihood'          }
    {'SmoothedStates'         }
    {'SmoothedStatesCov'      }
    {'SmoothedStateDisturb'   }
    {'SmoothedStateDisturbCov'}
    {'SmoothedObsInnov'       }
    {'SmoothedObsInnovCov'    }
    {'KalmanGain'             }
    {'DataUsed'               }

Convert Output to a table.

OutputTbl = struct2table(Output);
OutputTbl(1:10,1:4) % Display first ten rows of first four variables

ans=10×4 table
    LogLikelihood    SmoothedStates    SmoothedStatesCov    SmoothedStateDisturb
    _____________    ______________    _________________    ____________________

    {0x0 double}      {0x0 double}       {0x0 double}           {0x0 double}    
    {0x0 double}      {0x0 double}       {0x0 double}           {0x0 double}    
    {[  0.1827]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.0972]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.4472]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.2073]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.5167]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.2389]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[  0.5064]}      {2x1 double}       {2x2 double}           {2x1 double}    
    {[ -0.0105]}      {2x1 double}       {2x2 double}           {2x1 double}    

The first two rows of the table contain empty cells or zeros. These correspond to the observations required to initialize the diffuse Kalman filter. That is, SwitchTime is 2.

SwitchTime = 2;

Plot the smoothed states and their individual 95% Wald-type confidence intervals.

CI = nan(T,2,2);

for j = (SwitchTime + 1):T
    CovX = OutputTbl.SmoothedStatesCov{j};
    CI(j,:,1) = X(j,1) + 1.96*sqrt(CovX(1,1))*[-1 1];
    CI(j,:,2) = X(j,2) + 1.96*sqrt(CovX(2,2))*[-1 1];    
end

figure;
plot(1:T,X(:,1),'k',1:T,CI(:,:,1),'--r');
xlabel('Period');
ylabel('Smoothed states');
title('State 1 Estimates')
legend('Smoothed','95% Individual CIs');
grid on;

figure;
plot(1:T,X(:,2),'k',1:T,CI(:,:,2),'--r');
xlabel('Period');
ylabel('Smoothed states');
title('State 2 Estimates')
legend('Smoothed','95% Individual CIs');
grid on;