refine

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

Suppose further that the latent process is subject to additive measurement error. Consequently, the observation equation is

where $${\epsilon }_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 = x + randn(T,1);

Together, the latent process and observation equations compose a state-space model. Assume that the state is a stationary AR(1) process. Then the state-space model to estimate is

Specify the coefficient matrices. Use NaN values for unknown parameters.

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

Create the diffuse state-space model by passing the coefficient matrices to dssm and specifying that the state type is diffuse.

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

Mdl is a dssm model object. The software sets values for the initial state mean and variance to 0 and Inf. Verify that the model is specified correctly using the display in the Command Window.

Find a good set of starting parameters to use for estimation.

params0 = [1 1 1]; % Initial values chosen arbitrarily
Output = refine(Mdl,y,params0);

Output is a 1-by-5 structure array containing the recommended initial parameter values.

Choose the initial parameter values corresponding to the largest loglikelihood

logL = cell2mat({Output.LogLikelihood})';
[~,maxLogLIndx] = max(logL)

maxLogLIndx = 
3

refinedParams0 = Output(maxLogLIndx).Parameters

refinedParams0 = 1×3

    0.9781    0.8965    0.9336

Description = Output(maxLogLIndx).Description

Description = 
'Loose bound interior point'

The algorithm that yields the highest loglikelihood value is Quasi-Newton, which is the first struct in the structure array Output.

Estimate Mdl using refinedParams0, which is the vector of refined initial parameter values.

EstMdl = estimate(Mdl,y,refinedParams0,'lb',[-Inf,0,0]);

Method: Maximum likelihood (fmincon)
Effective Sample size:             99
Logarithmic  likelihood:     -179.018
Akaike   info criterion:      364.036
Bayesian info criterion:      371.851
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.97805       0.02947   33.18393   0      
 c(2) |  0.89651       0.18465    4.85529  0.00000 
 c(3) |  0.93355       0.15187    6.14707   0      
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | -3.95108       0.72269   -5.46719   0      

The AR model coefficient is within two standard errors of 1, which suggests that the state processes is a random walk.

Suppose that the relationship between the unemployment rate and the nominal gross national product (nGNP) is linear. Suppose further that the 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 data.

load Data_NelsonPlosser

Preprocess the data by taking the first difference of the unemployment rate and converting nGNP to a series of returns. Remove the observations corresponding to the sequence of NaN values at the beginning of the unemployment rate series.

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
y = DataTable.UR(~isNaN);
y = diff(y);
T = size(y,1);
Z = [ones(T,1) 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);

Mdl is a dssm model object.

Find a good set of starting parameters to use for estimation.

params0 = [150 1000]; % Initial values chosen arbitrarily
Beta0 = [1 -100];
Output = refine(Mdl,y,params0,'Predictors',Z,'Beta0',Beta0);

Output is a 1-by-5 structure array containing the recommended initial parameter values.

Choose the initial parameter values corresponding to the largest loglikelihood.

logL = cell2mat({Output.LogLikelihood})';
[~,maxLogLIndx] = max(logL)

maxLogLIndx = 
4

refinedParams0 = Output(maxLogLIndx).Parameters

refinedParams0 = 1×4

    0.2070    1.3229    1.3610  -24.8848

Description = Output(maxLogLIndx).Description

Description = 
'Starting value perturbation'

The algorithm that yields the highest loglikelihood value is Nelder-Mead simplex, which is the second struct in the structure array Output.

Estimate Mdl using the refined initial parameter values refinedParams0.

EstMdl = estimate(Mdl,y,refinedParams0(1:(end - 2)),'Predictors',Z,...
    'Beta0',refinedParams0((end - 1):end));

Method: Maximum likelihood (fminunc)
Effective Sample size:             60
Logarithmic  likelihood:     -101.924
Akaike   info criterion:      211.849
Bayesian info criterion:      220.292
           |      Coeff       Std Err    t Stat     Prob  
----------------------------------------------------------
 c(1)      |   0.20700       0.12330     1.67891  0.09317 
 c(2)      |   1.32287       0.08415    15.71964   0      
 y <- z(1) |   1.36101       0.23736     5.73388   0      
 y <- z(2) | -24.88484       1.78021   -13.97861   0      
           |                                              
           |    Final State   Std Dev     t Stat    Prob  
 x(1)      |   1.21611        0           Inf      0