refine

Class: ssm

Refine initial parameters to aid estimation of state-space models

Syntax

Description

example

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.

example

refine(Mdl,Y,params0,Name,Value) displays results of the routines with additional options specified by one or more Name,Value pair arguments.

example

Output = refine(___) returns a structure array (Output) containing a vector of refined, initial parameter values, the loglikelihood corresponding the initial parameter values, and the method the software used to obtain the values using any of the input arguments in the previous syntaxes.

Tips

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

  • The refined, initial parameter values returned by refine might 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 and its initial parameter values to [].

Input Arguments

expand all

Mdl — State-space modelssm model

State-space model containing unknown parameters, specified as an ssm model returned by ssm.

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 datacell vector of numeric vectors | matrix

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

  • If Mdl is time invariant with respect to the observation equation, then Y is a T-by-n matrix, where each row corresponds to a period and each column corresponds to a particular observation in the model. Therefore, T is the sample size and m 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. Each element of the cell vector corresponds to a period and contains an nt-dimensional vector of observations for that period. 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 you created Mdl implicitly by specifying a parameter-to-matrix mapping function, and the function has input arguments for the observed responses or predictors (see the ParamMap name-value pair argument of ssm), the function then establishes a link to observed responses and the predictor data. The link overrides the value of Y.

Data Types: double | cell

params0 — Unknown parameter initial valuesnumeric vector

Unknown parameter initial values for numerical 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, and Cov0.

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

Data Types: double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

'Beta0' — Initial values of regression coefficientsmatrix

Initial values of regression coefficients, specified as the comma-separated pair consisting of 'Beta0' and a matrix.

The number of rows in Beta0 must equal the number of columns in Y. Beta0 and Predictors must have the same number of columns.

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

Data Types: double

'Predictors' — Predictor variables in state-space model observation equation[] (default) | matrix

Predictor variables in the state-space model observation equation, specified as the comma-separated pair consisting of 'Predictors' and a matrix. The columns of Predictors correspond to individual predictor variables. Predictors must have T rows, where row t corresponds to the observed predictors at period t (Zt) in the expanded observation equation

ytZtβ=Cxt+Dut.

That is, the software deflates the observations using the regression component. β is the time-invariant vector of regression coefficients that the software estimates with all other parameters.

If there are n observations per period, then the software regresses all predictor series onto each observation. Subsequently, the software returns d-by-n matrix of fitted regression coefficient vectors for each observation series, where d is the number of columns of Predictors.

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

expand all

Output — Information about initial parameter valuesstructure 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.

FieldDescription
Description

Refinement algorithm.

Each element of Output corresponds to one of the following algorithms:

'Loose bound interior point'
'Nelder-Mead algorithm'
'Quasi-Newton'
'Starting value perturbation'
'Starting value shrinkage'
LoglikelihoodLoglikelihood corresponding to the initial parameter values.
ParametersVector 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 might improve.

Examples

expand all

Refine Parameters When Fitting a Time-Invariant, State-Space Model

Suppose that a latent process is a random walk. Subsequently, 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. Subsequently, the observation equation is

$$y_t = x_t + \varepsilon_t,$$

where $\varepsilon_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 = \phi x_{t-1} + \sigma_1u_t \\ y_t = x_t +
\sigma_2\varepsilon_t.\end{array}$$

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

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

Specify the state-space model using the coefficient matrices. Specify that the initial state distribution is stationary using the StateType name-value pair argument.

StateType = 0;
Mdl = ssm(A,B,C,D,'StateType',StateType);

Mdl is an ssm model. The software sets values for the initial state mean and variance. Verify that the model is correctly specified using the display in the Command Window.

Pass the observations to estimate to estimate the parameters. Set starting values for the parameters to params0 that are likely not their corresponding true values. Also, specify lower bound constraints of 0 for the standard deviations.

params0 = (1e-4)*ones(3,1);
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf,0,0]);
Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 3000 (the default value).

Warning: Numerical optimization does not return a convergence flag. Maximum
likelihood estimation is less likely to succeed if starting parameter values are
far from true values. Try different starting values, other optimization
algorithms and settings. Also try to reparameterize the model and rescale the
data. 

Method: Maximum likelihood (fmincon)
Sample size: 100
Logarithmic  likelihood:     -190.082
Akaike   info criterion:      386.164
Bayesian info criterion:       393.98
      |     Coeff       Std Err     t Stat      Prob  
------------------------------------------------------
 c(1) |  0.91794        0.04164     22.04345   0      
 c(2) |  1.35542        0.18920      7.16394   0      
 c(3) |  0.00416       24.44536      0.00017  0.99986 
      |                                               
      |   Final State    Std Dev      t Stat    Prob  
 x(1) | -3.37650        0.00416   -811.26251   0      

estimate failed to converge, and so the results are undesirable.

Refine params0 using refine.

Output = refine(Mdl,y,params0);
logL = cell2mat({Output.LogLikelihood})';
[~,maxLogLIndx] = max(logL)
refinedParams0 = Output(maxLogLIndx).Parameters
Description = Output(maxLogLIndx).Description
The likelihood is not well defined at the starting parameter values.

maxLogLIndx =

     3


refinedParams0 =

    0.9705
    0.8934
    0.9330


Description =

Loose bound interior point

The algorithm that yields the highest loglikelihood value is Loose bound interior point, which is the third 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)
Sample size: 100
Logarithmic  likelihood:     -181.379
Akaike   info criterion:      368.758
Bayesian info criterion:      376.574
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.97050       0.02863   33.90368   0      
 c(2) |  0.89343       0.18521    4.82401  0.00000 
 c(3) |  0.93303       0.15176    6.14806   0      
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | -3.93007       0.72066   -5.45343   0      

estimate converged, making the parameter estimates much more desirable. The AR model coefficient is within two standard errors of 1, which suggests that the state processes is a random walk.

Refine Estimation of a State-Space Model That Includes a Regression Component

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

$$\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{x_{1,t}}}\\
{{x_{2,t}}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
\phi &\theta \\
0&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_{1,t - 1}}}\\
{{x_{2,t - 1}}}
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
1&0\\
1&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{u_{1,t}}}\\
{{u_{2,t}}}
\end{array}} \right]\\
{y_t} - \beta {Z_t} = {x_{1,t}} + \sigma\varepsilon_t,
\end{array}$$

where:

  • $x_{1,t}$ is the change in the unemployment rate at time t.

  • $x_{2,t}$ is a dummy state for the MA(1) effect.

  • $y_{t}$ is the observed change in the unemployment rate being deflated by the growth rate of nGNP ( $Z_t$).

  • $u_{1,t}$ is the Gaussian series of state disturbances having mean 0 and standard deviation 1.

  • $\varepsilon_t$ is the Gaussian series of observation innovations having mean 0 and standard deviation $\sigma$.

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. 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);                          % The sample size
Z = [ones(T-1,1) diff(log(gnpn))];
y = diff(u);

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 NaN; 0 0];
B = [1 0;1 0];
C = [1 0];
D = NaN;

Specify the state-space model using ssm.

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

Estimate the model parameters. 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 0 1e-11];
Beta0 = [0 0];
EstMdl = estimate(Mdl,y,params0,'Predictors',Z,...
    'Beta0',Beta0,'lb',[-Inf,-Inf,0,-Inf,-Inf]);
Warning: Covariance matrix of estimators cannot be computed precisely due to
inversion difficulty. Check parameter identifiability. Also try different
starting values and other options to compute the covariance matrix. 

Method: Maximum likelihood (fmincon)
Sample size: 61
Logarithmic  likelihood:     -109.709
Akaike   info criterion:      229.417
Bayesian info criterion:      239.972
           |      Coeff        Std Err      t Stat     Prob  
-------------------------------------------------------------
 c(1)      |  -0.25172       0.27386       -0.91917  0.35801 
 c(2)      |   0.50899       0.23603        2.15650  0.03105 
 c(3)      |   0.00000       2.56261e+07    0.00000   1      
 y <- z(1) |   1.85671       0.11960       15.52385   0      
 y <- z(2) | -27.50973       1.01219      -27.17850   0      
           |                                                 
           |    Final State   Std Dev        t Stat    Prob  
 x(1)      |   0.84457        0              Inf      0      
 x(2)      |   0.88637        0              Inf      0      

The software could not estimate the parameter covariance matrix.

Refine the initial parameter values.

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)
refinedParams0 = Output(maxLogLIndx).Parameters
Description = Output(maxLogLIndx).Description
maxLogLIndx =

     1


refinedParams0 =

   -0.3410    1.0500    0.4859    1.3612  -24.4671


Description =

Quasi-Newton

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

Estimate Mdl using the refined initial parameter values refinedParams0.

pBeta = numel(Beta0);
EstMdl = estimate(Mdl,y,refinedParams0(1:(end - pBeta)),'Predictors',Z,...
    'Beta0',refinedParams0((end - pBeta + 1):end),...
    '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 

estimate returns reasonable parameter estimates and their corresponding standard errors.

Was this topic helpful?