refine

Refine initial parameters to aid diffuse state-space model estimation

Syntax

refine(Mdl,Y,params0)

refine(Mdl,Y,params0,Name,Value)

Output = refine(___)

Description

refine(Mdl,Y,params0) finds a set of initial parameter values to use when fitting the state-space model Mdl to the response data Y, using the crude set of initial parameter values params0. The software uses several routines, and displays the resulting loglikelihood and initial parameter values for each routine.

refine(Mdl,Y,params0,Name,Value) displays results of the routines with additional options specified by one or more Name,Value pair arguments. For example, you can include a linear regression component composed of predictors and an initial value for the coefficients.

Output = refine(___) returns a structure array (Output) containing a vector of refined, initial parameter values, the loglikelihood corresponding the initial parameter values, and the computation method yielding the values. You can use any of the input arguments in the previous syntaxes.

Input Arguments

Mdl — Diffuse state-space model

dssm model object

Diffuse state-space model containing unknown parameters, specified as a dssm model object returned by dssm.

• For explicitly created state-space models, the software estimates all NaN values in the coefficient matrices (Mdl.A, Mdl.B, Mdl.C, and Mdl.D) and the initial state means and covariance matrix (Mdl.Mean0 and Mdl.Cov0). For details on explicit and implicit model creation, see dssm.

• For implicitly created state-space models, you specify the model structure and the location of the unknown parameters using the parameter-to-matrix mapping function. Implicitly create a state-space model to estimate complex models, impose parameter constraints, and estimate initial states. The parameter-to-mapping function can also accommodate additional output arguments.

Note

Mdl does not store observed responses or predictor data. Supply the data wherever necessary using, the appropriate input and name-value pair arguments.

Y — Observed response data

numeric matrix | cell vector of numeric vectors

Observed response data to which Mdl is fit, specified as a numeric matrix or a cell vector of numeric vectors.

• If Mdl is time invariant with respect to the observation equation, then Y is a T-by-n matrix. Each row of the matrix corresponds to a period and each column corresponds to a particular observation in the model. Therefore, T is the sample size and n is the number of observations per period. The last row of Y contains the latest observations.

• If Mdl is time varying with respect to the observation equation, then Y is a T-by-1 cell vector. Y{t} contains an n_t-dimensional vector of observations for period t, where t = 1,...,T. The corresponding dimensions of the coefficient matrices in Mdl.C{t} and Mdl.D{t} must be consistent with the matrix in Y{t} for all periods. The last cell of Y contains the latest observations.

Suppose that you create Mdl implicitly by specifying a parameter-to-matrix mapping function, and the function has input arguments for the observed responses or predictors. Then, the mapping function establishes a link to observed responses and the predictor data in the MATLAB^® workspace, which overrides the value of Y.

NaN elements indicate missing observations. For details on how the Kalman filter accommodates missing observations, see Algorithms.

Data Types: double | cell

params0 — Initial values of unknown parameters

numeric vector

Initial values of unknown parameters for numeric maximum likelihood estimation, specified as a numeric vector.

The elements of params0 correspond to the unknown parameters in the state-space model matrices A, B, C, and D, and, optionally, the initial state mean Mean0 and covariance matrix Cov0.

• If you created Mdl explicitly (that is, by specifying the matrices without a parameter-to-matrix mapping function), then the software maps the elements of params to NaNs in the state-space model matrices and initial state values. The software searches for NaNs column-wise, following the order A, B, C, D, Mean0, Cov0.

• If you created Mdl implicitly (that is, by specifying the matrices with a parameter-to-matrix mapping function), then set initial parameter values for the state-space model matrices, initial state values, and state types within the parameter-to-matrix mapping function.

Data Types: double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Beta0 — Initial values of regression coefficients

numeric matrix

Initial values of regression coefficients, specified as the comma-separated pair consisting of 'Beta0' and a d-by-n numeric matrix. d is the number of predictor variables (see Predictors) and n is the number of observed response series (see Y).

By default, Beta0 is the ordinary least-squares estimate of Y onto Predictors.

Data Types: double

Predictors — Predictor data

[] (default) | numeric matrix

Predictor data for the regression component in the observation equation, specified as the comma-separated pair consisting of 'Predictors' and a T-by-d numeric matrix. T is the number of periods and d is the number of predictor variables. Row t corresponds to the observed predictors at period t (Z_t) in the expanded observation equation

$$y_t-Z_t\beta =C{x}_t+D{u}_t.$$

In other words, the predictor series serve as observation deflators. β is a d-by-n time-invariant matrix of regression coefficients that the software estimates with all other parameters.

• For n observations per period, the software regresses all predictor series onto each observation.

• If you specify Predictors, then Mdl must be time invariant. Otherwise, the software returns an error.

• By default, the software excludes a regression component from the state-space model.

Data Types: double

Output Arguments

Output — Information about initial parameter values

structure array

Information about the initial parameter values, returned as a 1-by-5 structure array. The software uses five algorithms to find initial parameter values, and each element of Output corresponds to an algorithm.

This table describes the fields of Output.

Field	Description
Description	Refinement algorithm. Each element of Output corresponds to one of these algorithms: 'Loose bound interior point' 'Nelder-Mead algorithm' 'Quasi-Newton' 'Starting value perturbation' 'Starting value shrinkage'
Loglikelihood	Loglikelihood corresponding to the initial parameter values.
Parameters	Vector of refined initial parameter values. The order of the parameters is the same as the order in params0. If you pass these initial values to estimate, then the estimation results can improve.

Examples

Refine Parameters When Fitting Time-Invariant Diffuse State-Space Model

Suppose that a latent process is a random walk. Consequently, the state equation is

$$x_t=x_{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; 
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

$$y_t=x_t+{\epsilon }_t,$$

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

$$\begin{array}{c}{x}_t=\varphi x_{t-1}+{\sigma }_1{u}_t\\ y_t=x_t+{\sigma }_2{\epsilon }_t.\end{array}$$

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.

Refine Diffuse State-Space Model Estimation Including Regression Component

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

$$\begin{array}{l}{x}_t=\varphi x_{t-1}+\sigma u_t\\ y_t-\beta Z_t=x_t,\end{array}$$

where:

• $$x_t$$ is the unemployment rate at time t.

• $$y_t$$ is the observed unemployment rate being deflated by the log of nGNP ($$Z_t$$).

• $$u_t$$ is the Gaussian series of state disturbances having mean 0 and unknown standard deviation $$\sigma$$.

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      

Tips

• Likelihood surfaces of state-space models can be complicated, for example, they can contain multiple local maxima. If estimate fails to converge, or converges to an unsatisfactory solution, then refine can find a better set of initial parameter values to pass to estimate.

• The refined initial parameter values returned by refine can appear similar to each other and to params0. Choose a set yielding estimates that make economic sense and correspond to relatively large loglikelihood values.

• If a refinement attempt fails, then the software displays errors and sets the corresponding loglikelihood to -Inf. It also sets its initial parameter values to [].

Algorithms

The Kalman filter accommodates missing data by not updating filtered state estimates corresponding to missing observations. In other words, suppose that your data has a missing observation at period t. Then, the state forecast for period t, based on the previous t – 1 observations, is equivalent to the filtered state for period t.

Version History

Introduced in R2015b