Vector Autoregressive with External inputs (VARX)

In some fields, the use of VARX is widespread. VARX models are special cases of general linear models on the form

\[\begin{aligned} x^+ &= Ax + Bu + Ke\\ y &= Cx + Du + e \end{aligned}\]

which we show here by simulating the VARX model

\[y_t = A_1 y_{t-1} + A_2 y_{t-2} + B_1 u_{t-1}\]

and estimating a regular statespace model.

using ControlSystemIdentification, Plots

A1 = 0.3randn(2,2)
A2 = 0.3randn(2,2)
B1 = randn(2,2)

N = 300
Y = [randn(2), randn(2)]
u = randn(2, N)

for i = 3:N
    yi = A1*Y[i-1] + A2*Y[i-2] + B1*u[:,i-1]
    push!(Y, yi)
end

y = hcat(Y...)

d = iddata(y, u, 1)

InputOutput data of length 300, 2 outputs, 2 inputs, Ts = 1

We now estimate two models, one using subspace-based identification (subspaceid) and one using the prediction-error method (newpem). The VARX model we estimated had a state of order 4, two lags of $y$, each of which is a vector of length 2, we thus estimate models of order 4 below.

model1 = subspaceid(d, 4, zeroD=true)
model2, x0 = newpem(d, 4)

plot(
    simplot(d, model1, title="Simulation subspaceid", layout=2),
    simplot(d, model2, x0, title="Simulation PEM", layout=2),
    titlefontsize=10,
)

The simulation indicates that the fit is close to 100%, i.e., the general linear model fit the VARX model perfectly, it does not however have exactly the same structure as the original VARX model.

The estimated models on statespace form can be converted to MIMO transfer functions (polynomial models) by calling tf:

using ControlSystemsBase: tf
tf(model2.sys)

TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}}
Input 1 to output 1
 0.038289788045511586z^3 + 0.0019646642624622235z^2 - 0.03595727331913012z - 3.4679203952947546e-13
-----------------------------------------------------------------------------------------------------
1.0z^4 + 0.39655401871518237z^3 - 0.31834441483598674z^2 - 0.13111984528997145z - 0.11677055287897602

Input 1 to output 2
   -0.39153150668015124z^3 - 0.04590879957615146z^2 + 0.21539652003408521z + 6.224604165439018e-13
-----------------------------------------------------------------------------------------------------
1.0z^4 + 0.39655401871518237z^3 - 0.31834441483598674z^2 - 0.13111984528997145z - 0.11677055287897602

Input 2 to output 1
   -1.405109762210532z^3 - 0.39176446503107754z^2 - 0.24202080721627428z + 2.7670921110001245e-13
-----------------------------------------------------------------------------------------------------
1.0z^4 + 0.39655401871518237z^3 - 0.31834441483598674z^2 - 0.13111984528997145z - 0.11677055287897602

Input 2 to output 2
   0.10183534775337588z^3 + 0.028182914569904627z^2 - 0.3241332547868683z - 6.252220963176569e-13
-----------------------------------------------------------------------------------------------------
1.0z^4 + 0.39655401871518237z^3 - 0.31834441483598674z^2 - 0.13111984528997145z - 0.11677055287897602

Sample Time: 1.0 (seconds)
Discrete-time transfer function model

Summary

The statespace model

\[\begin{aligned} x^+ &= Ax + Bu + Ke\\ y &= Cx + Du + e \end{aligned}\]

is very general and subsumes a lot of other, special-structure models, like the VARX model.