Estimate parameters of ARMAX model using time-domain data

returns an `sys`

= armax(`data`

,```
[na
nb nc nk]
```

)`idpoly`

model, `sys`

,
with estimated parameters and covariance (parameter uncertainties).
Estimates the parameters using the prediction-error method and specified
polynomial orders.

returns
an `sys`

= armax(`data`

,```
[na
nb nc nk]
```

,`Name,Value`

)`idpoly`

model, `sys`

, with
additional options specified by one or more `Name,Value`

pair
arguments.

`data`

— Time-domain estimation data`iddata`

objectTime-domain estimation data, specified as an `iddata`

object.
You cannot use frequency-domain data for estimating ARMAX models.

`[na nb nc nk]`

— Polynomial orders1-by-4 vector of positive integers | 1-by-4 vector of matricesPolynomial orders of an ARMAX Model,
specified as a 1-by-4 vector, `[na nb nc nk]`

.

For a model with *Ny* outputs and *Nu* inputs:

`na`

is the order of the polynomial*A*(*q*), specified as an*Ny*-by-*Ny*matrix of nonnegative integers.`nb`

is the order of the polynomial*B*(*q*) + 1, specified as an*Ny*-by-*Nu*matrix of nonnegative integers.`nc`

is the order of the polynomial*C*(*q*), specified as a column vector of nonnegative integers of length*Ny*.`nk`

is the input-output delay expressed as fixed leading zeros of the*B*polynomial.Specify

`nk`

as an*Ny*-by-*Nu*matrix of nonnegative integers.

`init_sys`

— System for configuring initial parametrizationdiscrete-time linear modelSystem for configuring initial parametrization of `sys`

,
specified as a discrete-time linear model. You obtain `init_sys`

by
either performing an estimation using measured data or by direct construction
using commands such as `idpoly`

and `idss`

.

If `init_sys`

is an ARMAX model, `armax`

uses
the parameter values of `init_sys`

as the initial
guess for estimation. To configure initial guesses and constraints
for *A*(*q*), *B*(*q*),
and *C*(*q*), use the `Structure`

property
of `init_sys`

. For example:

To specify an initial guess for the

*A*(*q*) term of`init_sys`

, set`init_sys.Structure.a.Value`

as the initial guess.To specify constraints for the

*B*(*q*) term of`init_sys`

:set

`init_sys.Structure.b.Minimum`

to the minimum*B*(*q*) coefficient values.set

`init_sys.Structure.b.Maximum`

to the maximum*B*(*q*) coefficient values.set

`init_sys.Structure.b.Free`

to indicate which*B*(*q*) coefficients are free for estimation.

If `init_sys`

is not a polynomial model of
ARMAX structure, the software first converts `init_sys`

to
an ARMAX model. `armax`

uses the parameters of
the resulting model as the initial guess for estimating `sys`

.

If `opt`

is not specified, and `init_sys`

was
obtained by estimation, then the estimation options from `init_sys.Report.OptionsUsed`

are
used.

`opt`

— Estimation options`armaxOptions`

option setEstimation options for ARMAX model identification, specified
as an `armaxOptions`

option set.

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`0`

(default) | scalar | vector of positive integersInput delays, specified as the comma-separated pair consisting
of `'InputDelay'`

and one of the following:

`Nu`

-by-1 vector, where`Nu`

is the number of inputs — Each entry is a numerical value representing the input delay for the corresponding input channel. Specify input delays as integer multiples of the sample time`Ts`

. For example,`InputDelay = 3`

means a delay of three sampling periods.Scalar value — The same delay is applied to all input channels.

`'ioDelay'`

— Transport delays`0`

(default) | scalar | matrixTransport delays for each input/output pair, specified as the
comma-separated pair consisting of `'ioDelay'`

and
one of the following:

`Ny`

-by-`Nu`

matrix, where`Ny`

is the number of outputs and`Nu`

is the number of inputs — Each entry is an integer value representing the transport delay for the corresponding input/output pair. Specify transport delays as integers denoting delay of a multiple of the sample time`Ts`

.Scalar value — The same delay is applied to all input/output pairs.

`'ioDelay'`

is useful as a replacement for
the `nk`

order. You can factor out `max(nk-1,0)`

lags
as the `ioDelay`

value. For `nk`

>1, `armax(na,nb,nk)`

is
equivalent to `armax(na,nb,1,'ioDelay',nk-1)`

.

`'IntegrateNoise'`

— Addition of integrators in noise channel`false(Ny,1)`

, where Addition of integrators in noise channel, specified as the comma-separated
pair consisting of `‘IntegrateNoise'`

and
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(q)y(t)=B(q)u(t-nk)+\frac{C(q)}{1-{q}^{-1}}e(t)$$

Where, $$\frac{1}{1-{q}^{-1}}$$ is the integrator in the noise
channel, *e*(*t*).

Use `IntegrateNoise`

to create ARIMA or ARIMAX
models.

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

`sys`

— ARMAX model`idpoly`

objectARMAX model that fits the given estimation data, returned as
a discrete-time `idpoly`

object.
This model is created using the specified model orders, delays, and
estimation options.

Information about the estimation results and options used is
stored in the `Report`

property of the model. `Report`

has
the following fields:

Report Field | Description | ||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

`Status` | Summary of the model status, which indicates whether the model was created by construction or obtained by estimation. | ||||||||||||||||||

`Method` | Estimation command used. | ||||||||||||||||||

`InitialCondition` | Handling of initial conditions during model estimation, returned as a string with one of the following values: `'zero'` — The initial conditions were set to zero.`'estimate'` — The initial conditions were treated as independent estimation parameters.`'backcast'` — The initial conditions were estimated using the best least squares fit.
This field is especially useful to view
how the initial conditions were handled when the | ||||||||||||||||||

`Fit` | Quantitative assessment of the estimation, returned as a structure. See Loss Function and Model Quality Metrics for more information on these quality metrics. The structure has the following fields:
| ||||||||||||||||||

`Parameters` | Estimated values of model parameters. | ||||||||||||||||||

`OptionsUsed` | Option set used for estimation. If no custom options
were configured, this is a set of default options. See | ||||||||||||||||||

`RandState` | State of the random number stream at the start of estimation.
Empty, | ||||||||||||||||||

`DataUsed` | Attributes of the data used for estimation, returned as a structure with the following fields:
| ||||||||||||||||||

`Termination` | Termination conditions for the iterative search used for prediction error minimization. Structure with the following fields:
For estimation methods that
do not require numerical search optimization, the |

For more information on using `Report`

, see Estimation Report.

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.

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

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

Load the estimation data.

load iddata2 z2

Estimate a state-space model of order 3 from the estimation data.

sys0 = n4sid(z2,3);

Estimate an ARMAX model using the previously estimated state-space model as an initial guess.

sys = armax(z2,sys0);

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

`armax`

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

to
estimate an Output-Error (OE) model.

The ARMAX model structure is:

$$\begin{array}{l}y(t)+{a}_{1}y(t-1)+\dots +{a}_{{n}_{a}}y(t-{n}_{a})=\\ \text{}{b}_{1}u(t-{n}_{k})+\dots +{b}_{{n}_{b}}u(t-{n}_{k}-{n}_{b}+1)+\\ \text{}{c}_{1}e(t-1)+\dots +{c}_{{n}_{c}}e(t-{n}_{c})+e(t)\end{array}$$

A more compact way to write the difference equation is:

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

where,

$$y(t)$$— 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(t-1)\dots y(t-{n}_{a})$$ — Previous outputs on which the current output depends.

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

$$e(t-1)\dots e(t-{n}_{c})$$ — 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(q)=1+{a}_{1}{q}^{-1}+\dots +{a}_{{n}_{a}}{q}^{-{n}_{a}}$$

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

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

If `data`

is a time series that has no input
channels and one output channel, then `armax`

calculates
an ARMA model for the time series

$$A(q)y(t)=C(q)e(t)$$

In this case,

orders = [na nc]

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

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

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

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

Use the

`IntegrateNoise`

property to add integrators to the noise source.

An iterative search algorithm minimizes a robustified quadratic prediction error criterion. The iterations are terminated when any of the following is true:

Maximum number of iterations is reached.

Expected improvement is less than the specified tolerance.

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(q)$$ and $$F(q)$$ (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.

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

`aic`

| `armaxOptions`

| `arx`

| `bj`

| `forecast`

| `fpe`

| `goodnessofFit`

| `iddata`

| `idfrd`

| `idpoly `

| `oe`

| `polyest`

| `ssest`

| `tfest`

Was this topic helpful?