# armax

Estimate parameters of ARMAX model using time-domain data

## Syntax

```sys = armax(data,[na nb nc nk])sys = armax(data,[na nb nc nk],Name,Value)sys = armax(data,init_sys)sys = armax(data,___,opt)```

## Description

 Note:   `armax` supports only time-domain data. For frequency-domain data, use `oe`.

```sys = armax(data,[na nb nc nk])``` returns an `idpoly` model, `sys`, with estimated parameters and covariance (parameter uncertainties). Estimates the parameters using the prediction-error method and specified polynomial orders.

```sys = armax(data,[na nb nc nk],Name,Value)``` returns an `idpoly` model, `sys`, with additional options specified by one or more `Name,Value` pair arguments.

`sys = armax(data,init_sys)` estimates a polynomial model using the ARMAX structure polynomial model `init_sys` to configure the initial parameterization.

`sys = armax(data,___,opt)` specifies estimation options using the option set `opt`.

## Input Arguments

 `data` Estimation data. Specify `data` as an `iddata` object containing the time-domain input-output data. You cannot use frequency-domain data for estimating ARMAX models. `[na nb nc nk]` Polynomial orders. `[na nb nc nk]` define the polynomial orders of an ARMAX Model. `na` — Order of the polynomial A(q). Specify `na` as an Ny-by-Ny matrix of nonnegative integers. Ny is the number of outputs.`nb` — Order of the polynomial B(q) + 1.`nb` is an Ny-by-Nu matrix of nonnegative integers. Ny is the number of outputs and Nu is the number of inputs.`nc` — Order of the polynomial C(q).`nc` is a column vector of nonnegative integers of length Ny. Ny is the number of outputs.`nk` — Input-output delay expressed as fixed leading zeros of the B polynomial. Specify `nk` as an Ny-by-Nu matrix of nonnegative integers. Ny is the number of outputs and Nu is the number of inputs. `init_sys` Linear polynomial model that configures the initial parameterization of `sys`. `init_sys` must be an ARMAX model. You may obtain `init_sys` by either performing an estimation using measured data, or by direct construction. Use the `Structure` property of `init_sys` to configure initial guesses and constraints for A(q), B(q), and C(q). To specify an initial guess for, say, the A(q) term of `init_sys`, set `init_sys.Structure.a.Value` as the initial guess. To specify constraints for, say, the B(q) term of `init_sys`: set `init_sys.Structure.b.Minimum` to the minimum B(q) coefficient valuesset `init_sys.Structure.b.Maximum` to the maximum B(q) coefficient valuesset `init_sys.Structure.b.Free` to indicate which B(q) coefficients are free for estimation You can similarly specify the initial guess and constraints for the other polynomials. If `opt` is not specified, and `init_sys` was created by estimation, then the estimation options from `init_sys.Report.OptionsUsed` are used. `opt` Estimation options. `opt` is an options set that specifies estimation options, including: estimation objectivehandling of initial conditionsnumerical search method to be used in estimation Use `armaxOptions` to create the options set.

### 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`.

 `'InputDelay'` Input delays. `InputDelay` is a numeric vector specifying a time delay for each input channel. Specify input delays in integer multiples of the sample time `Ts`. For example, `InputDelay = 3` means a delay of three sampling periods. For a system with `Nu` inputs, set `InputDelay` to an `Nu`-by-1 vector, where each entry is a numerical value representing the input delay for the corresponding input channel. You can also set `InputDelay` to a scalar value to apply the same delay to all channels. Default: 0 for all input channels `'ioDelay'` Transport delays. `ioDelay` is a numeric array specifying a separate transport delay for each input/output pair. Specify transport delays as integers denoting delay of a multiple of the sample time `Ts`. For a MIMO system with `Ny` outputs and `Nu` inputs, set `ioDelay` to a `Ny`-by-`Nu` array, where each entry is a numerical value representing the transport delay for the corresponding input/output pair. You can also set `ioDelay` to a scalar value to apply the same delay to all input/output pairs. Useful as a replacement for the `nk` order, you can factor out `max(nk-1,0)` lags as the `ioDelay` value. Default: 0 for all input/output pairs `'IntegrateNoise'` Logical vector specifying integrators in the noise channel. `IntegrateNoise` is a logical vector of length Ny, where Ny is the number of outputs. Setting `IntegrateNoise` to `true` for a particular output results in the model: $A\left(q\right)y\left(t\right)=B\left(q\right)u\left(t-nk\right)+\frac{C\left(q\right)}{1-{q}^{-1}}e\left(t\right)$ Where, $\frac{1}{1-{q}^{-1}}$ is the integrator in the noise channel, e(t). Use `IntegrateNoise` to create an ARIMA model. For example, ```load iddata9 z9; z9.y = cumsum(z9.y); %integrated data sys = armax(z9,[4 1],'IntegrateNoise',true); compare(z9,sys,10) %10-step ahead prediction``` Default: `false(Ny,1)` (Ny is the number of outputs.)

## Output Arguments

 `sys` Identified ARMAX structure polynomial model. `sys` is a discrete-time `idpoly` model, which encapsulates the estimated A, B and C polynomials and the parameter covariance information.

## Examples

collapse all

### Estimate ARMAX Model Using Regularization

Estimate a regularized ARMAX model by converting a regularized ARX model.

Load data.

```load regularizationExampleData.mat m0simdata; ```

Estimate an unregularized ARMAX model of order 15.

```m1 = armax(m0simdata(1:150), [30 30 30 1]); ```

Estimate a regularized ARMAX model by determining Lambda value by trial and error.

```opt = armaxOptions; opt.Regularization.Lambda = 1; m2 = armax(m0simdata(1:150), [30 30 30 1], opt); ```

Obtain a lower-order ARMAX model by converting a regularized ARX model followed by order reduction.

```opt1 = arxOptions; [L,R] = arxRegul(m0simdata(1:150), [30 30 1]); opt1.Regularization.Lambda = L; opt1.Regularization.R = R; m0 = arx(m0simdata(1:150),[30 30 1],opt1); mr = idpoly(balred(idss(m0),7)); ```

Compare the model outputs against data.

```compare(m0simdata(150:end), m1, m2, mr, compareOptions('InitialCondition','z')); ```

### Specify Estimation Options

Estimate an ARMAX model from measured data and specify the estimation options.

Estimate an ARMAX model with simulation focus, using `'lm'` as the search method and maximum number of search iterations set to 10.

```load twotankdata; z = iddata(y,u,0.2); opt = armaxOptions; opt.Focus = 'simulation'; opt.SearchMethod = 'lm'; opt.SearchOption.MaxIter = 10; opt.Display = 'on'; sys = armax(z, [2 2 2 1], opt); ```

The termination conditions for measured component of the model shown in the progress viewer is that the maximum number of iterations were reached.

To improve results, re-estimate the model using a greater value for `MaxIter`, or continue iterations on the previously estimated model as follows:

```sys2 = armax(z, sys); compare(z, sys, sys2) ```

where `sys2` refines the parameters of `sys` to improve the fit to data.

### Estimate an ARIMA Model

Estimate a 4th order ARIMA model for univariate time-series data.

```load iddata9; z9.y = cumsum(z9.y); % integrated data model = armax(z9, [4 1], 'IntegrateNoise', true); compare(z9, model, 10) % 10-step ahead prediction ```

### Estimate ARMAX Models Iteratively

Estimate ARMAX models of varying orders iteratively from measured data.

Estimate ARMAX models of orders varying between 1 and 4 for dryer data.

```load dryer2; z = iddata(y2,u2,0.08,'Tstart',0); na = 2:4; nc = 1:2; nk = 0:2; models = cell(1,18); ct = 1; for i = 1:3 na_ = na(i); nb_ = na_; for j = 1:2 nc_ = nc(j); for k = 1:3 nk_ = nk(k); models{ct} = armax(z, [na_, nb_, nc_, nk_]); ct = ct+1; end end end ```

Stack the estimated models and compare their simulated responses to estimation data `z`.

```models = stack(1,models{:}); compare(z,models) ```

## Alternatives

`armax` does not support continuous-time model estimation. Use `tfest` to estimate a continuous-time transfer function model, or `ssest` to estimate a continuous-time state-space model.

## More About

collapse all

### ARMAX Model

The ARMAX model structure is

A more compact way to write the difference equation is

$A\left(q\right)y\left(t\right)=B\left(q\right)u\left(t-{n}_{k}\right)+C\left(q\right)e\left(t\right)$

where

• $y\left(t\right)$— Output at time $t$.

• ${n}_{a}$ — Number of poles.

• ${n}_{b}$ — Number of zeroes plus 1.

• ${n}_{c}$ — Number of C coefficients.

• ${n}_{k}$ — Number of input samples that occur before the input affects the output, also called the dead time in the system.

• $y\left(t-1\right)\dots y\left(t-{n}_{a}\right)$ — Previous outputs on which the current output depends.

• $u\left(t-{n}_{k}\right)\dots u\left(t-{n}_{k}-{n}_{b}+1\right)$ — Previous and delayed inputs on which the current output depends.

• $e\left(t-1\right)\dots e\left(t-{n}_{c}\right)$ — White-noise disturbance value.

The parameters `na`, `nb`, and `nc` are the orders of the ARMAX model, and `nk` is the delay. q is the delay operator. Specifically,

$A\left(q\right)=1+{a}_{1}{q}^{-1}+\dots +{a}_{{n}_{a}}{q}^{-{n}_{a}}$

$B\left(q\right)={b}_{1}+{b}_{2}{q}^{-1}+\dots +{b}_{{n}_{b}}{q}^{-{n}_{b}+1}$

$C\left(q\right)=1+{c}_{1}{q}^{-1}+\dots +{c}_{{n}_{c}}{q}^{-{n}_{c}}$

If `data` is a time series, which has no input channels and one output channel, then `armax` calculates an ARMA model for the time series

$A\left(q\right)y\left(t\right)=e\left(t\right)$

In this case

```orders = [na nc] ```

### ARIMAX Model

An ARIMAX model structure is similar to ARMAX, except that it contains an integrator in the noise source e(t):

$A\left(q\right)y\left(t\right)=B\left(q\right)u\left(t-nk\right)+\frac{1}{\left(1-{q}^{-1}\right)}e\left(t\right)$

If there are no inputs, this reduces to an ARIMA model:

$A\left(q\right)y\left(t\right)=\frac{1}{\left(1-{q}^{-1}\right)}e\left(t\right)$

### Tips

• Use the `IntegrateNoise` property to add integrators to the noise source.

### Algorithms

An iterative search algorithm minimizes a robustified quadratic prediction error criterion. The iterations are terminated either when the maximum number of iterations is reached, or when the expected improvement is less than the specified tolerance, or when a lower value of the criterion cannot be found. You can get information about the stopping criteria using `sys.Report.Termination`.

Use the `armaxOptions` option set to create and configure options affecting the estimation results. In particular, set the search algorithm attributes, such as `MaxIter` and `Tolerance`, using the `'SearchOption'` property.

When you do not specify initial parameter values for the iterative search as an initial model, they are constructed in a special four-stage LS-IV algorithm.

The cutoff value for the robustification is based on the `Advanced.ErrorThreshold` estimation option and on the estimated standard deviation of the residuals from the initial parameter estimate. It is not recalculated during the minimization. By default, no robustification is performed; the default value of `ErrorThreshold` option is 0.

To ensure that only models corresponding to stable predictors are tested, the algorithm performs a stability test of the predictor. Generally, both $C\left(q\right)$ and $F\left(q\right)$ (if applicable) must have all zeros inside the unit circle.

Minimization information is displayed on the screen when the estimation option `'Display'` is `'On'` or `'Full'`. With `'Display'` `='Full'`, both the current and the previous parameter estimates are displayed in column-vector form, listing parameters in alphabetical order. Also, the values of the criterion function (cost) are given and the Gauss-Newton vector and its norm are also displayed. With ```'Display' = 'On'``` only the criterion values are displayed.

## References

Ljung, L. System Identification: Theory for the User, Upper Saddle River, NJ, Prentice-Hal PTR, 1999. See chapter about computing the estimate.

## See Also

Was this topic helpful?

Get trial now