IEEE TRANSACTIONS ON AUTOMATIC CONTROL

1

Parametric Identification Using
Weighted Null-Space Fitting

arXiv:1708.03946v3 [] 26 Mar 2018

Miguel Galrinho, Cristian R. Rojas, Member, IEEE, and Håkan Hjalmarsson, Fellow, IEEE

the aforementioned issues. We will not attempt to fully review
this vast field, but below we highlight some of the milestones.
Instrumental variable (IV) methods [3] allow consistency
to be obtained in a large variety of settings without the
issue of non-convexity. Asymptotic efficiency can be obtained
for some problems using iterative algorithms [4,5]. However,
IV methods cannot achieve the CR bound with closed-loop
data [6].
Realization-based methods [7], which later evolved into
subspace methods [8], are non-iterative and thus attractive
for their computational efficiency. The bias issue for closedloop data has been overcome by more recent algorithms [9]–
[12]. However, structural information is difficult to incorporate,
and—even if a complete analysis is still unavailable (significant contributions have been provided [12]–[14])—subspace
methods are in general not believed to be as accurate as PEM.
Index Terms—System identification, least squares.
Some methods are based on fixing some parameters in
certain places of the cost function but not others to obtain
a quadratic cost function, so that the estimate can be obtained
I. I NTRODUCTION
by (weighted) least squares. Then, the fixed coefficients are
For parametric system identification, the prediction error
replaced by an estimate from the previous iteration in the
method (PEM) is the reference in the field. With open-loop weighting or in a filtering step. This leads to iterative methdata, consistency is guaranteed if the model can describe the
ods, which date back to [15]. Some of these methods have
system dynamics, irrespective of the used noise model. For
been denoted iterative quadratic maximum likelihood (IQML),
Gaussian noise1 and with a noise model able to describe the originally developed for filter design [16,17] and later applied
noise spectrum, PEM with a quadratic cost function correto dynamical systems [18]–[20]. Another classical example is
sponds to maximum likelihood (ML), and is asymptotically
the Steiglitz-McBride method [21] for estimating output-error
efficient with respect to the used model structure [1], meaning models, which is equivalent to IQML for an impulse-input
that the covariance of the estimate asymptotically achieves the
case [22]. In the identification field, weightings or filterings
Cramér-Rao (CR) lower bound: the best possible covariance have not been determined by statistical considerations. In
for consistent estimators.
this perspective, the result in [23], showing that the SteiglitzThere are two issues that may hinder successful application McBride method is not asymptotically efficient, is not surprisof PEM. The first—and most critical—is the risk of converging ing.
to a non-global minimum of the cost function, which is
Another approach is to estimate, in an intermediate step, a
in general not convex. Thus, PEM requires local non-linear more flexible model, followed by a model reduction step to
optimization algorithms and good initialization points. The recover a model with the desired structure. The motivation for
second issue concerns closed-loop data. In this case, PEM this procedure is that, in some cases, each step corresponds
is biased unless the noise model is flexible enough. For to a convex optimization problem or a numerically reliable
asymptotic efficiency, the noise model must be of correct order procedure. To guarantee asymptotic efficiency, it is important
and estimated simultaneously with the dynamic model.
that the intermediate model is a sufficient statistic and the
During the half decade since the publication of [2], alterna- model reduction step is performed in a statistically sound way.
tives to PEM/ML have appeared, addressing one or both of Indirect PEM [24] formalizes the requirements starting with
an over-parametrized model of fixed order and uses ML in the
Automatic Control Lab and ACCESS Linnaeus Center, School of Electrical
Engineering, KTH Royal Institute of Technology, SE-100 44 Stockholm, model reduction step. The latter step corresponds in general
Sweden. (e-mail: {galrinho, crro, hjalmars}@kth.se.)
to a weighted non-linear least-squares problem.
This work was supported by the Swedish Research Council under contracts
It has also been recognized that the intermediate model does
2015-05285 and 2016-06079.
not need to capture the true system perfectly, but only with
1 When maximum likelihood and asymptotic efficiency are discussed in the
sufficient accuracy. Subspace algorithms can be interpreted in
following, the standard assumption is that the noise is Gaussian.

Abstract—In identification of dynamical systems, the prediction error method using a quadratic cost function provides
asymptotically efficient estimates under Gaussian noise and
additional mild assumptions, but in general it requires solving a
non-convex optimization problem. An alternative class of methods
uses a non-parametric model as intermediate step to obtain the
model of interest. Weighted null-space fitting (WNSF) belongs to
this class. It is a weighted least-squares method consisting of three
steps. In the first step, a high-order ARX model is estimated. In a
second least-squares step, this high-order estimate is reduced to
a parametric estimate. In the third step, weighted least squares
is used to reduce the variance of the estimates. The method is
flexible in parametrization and suitable for both open- and closedloop data. In this paper, we show that WNSF provides estimates
with the same asymptotic properties as PEM with a quadratic
cost function when the model orders are chosen according to the
true system. Also, simulation studies indicate that WNSF may
be competitive with state-of-the-art methods.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

this way: for example, SSARX [11] estimates an ARX model
followed by a singular-value-decomposition (SVD) step and
least-squares estimation. For spectral estimation, the THREElike approach [25] is also a two-step procedure that first
obtains a non-parametric spectral estimate and then reduces
it to a parametric estimate that in some sense is closest to
the non-parametric one, and whose optimization function is
convex [26].
In the field of time-series estimation, methods based on an
intermediate high-order time series have also been suggested
as alternative to ML, whose properties were studied in [27].
Durbin’s first method [28] for auto-regressive moving-average
(ARMA) time series uses an intermediate high-order AR time
series to simulate the innovations sequence, which allows
obtaining the ARMA parameters by least squares. This method
does not have the same asymptotic properties as ML, unlike
Durbin’s second method [28]. The latter is an extension of [29],
which had been proposed for MA time series, whose parameters can be estimated from the high-order AR estimates using
least squares, with an accuracy that can be made arbitrarily
close to the Cramér-Rao bound by increasing the AR-model
order. When applied to ARMA time series, the idea to achieve
efficiency is to iterate between estimating the AR and MA
polynomials using this procedure, initialized with Durbin’s
first method. Another way to achieve efficiency from Durbin’s
first method as starting point was proposed in [30] by using an
additional filtering step with the MA estimates from Durbin’s
first method, and then re-estimating the ARMA-parameters.
The asymptotic properties of these methods have been
analyzed by considering the high order tending to infinity, but
“small” compared to the sample size. A preferable analysis
should handle the relation between the high order and the
sample size formally, as done in [31] to prove consistency
of the method in [30], where the high order is assumed to
tend to infinity as function of the sample size at a particular
rate. This class of methods has become popular for vector
ARMA time series, with several available algorithms using
different procedures for obtaining the asymptotic efficient
model parameter estimates from the estimated innovations
(e.g., [32]–[34], with further information in references in these
papers). Despite sharing the same asymptotic properties, which
have been analyzed with the high order as a function of the
sample size, these algorithms may have different computational requirements and finite sample properties.
For identification of dynamical systems, instead of using
the high-order model to estimate the innovations, it has been
suggested that identification of the model of interest can be
done by applying asymptotic ML directly to the high-order
model [35]. The ASYM method [36] is an instantiation of this
approach. Because an ARX-model estimate and its covariance
constitute a sufficient statistic as the model order grows,
this approach can produce asymptotically efficient estimates.
However, the plant and noise models are estimated separately,
preventing asymptotic efficiency for closed-loop data. Also,
although such model reduction procedures may have numerical
advantages over direct application of PEM [36], this approach
still requires local non-linear optimization techniques. The
Box-Jenkins Steiglitz-McBride (BJSM) method [37] instead

2

uses the Steiglitz-McBride method in the model reduction step,
resulting in asymptotically efficient estimates of the plant in
open loop. Two drawbacks of BJSM are that the number of
iterations is required to tend to infinity (as for the SteiglitzMcBride method) and that, similarly to [35] and [36], the CR
bound cannot be attained in closed loop. The Model Order
Reduction Steiglitz-McBride (MORSM) method solves the
first drawback of BJSM, but not the second [38].
In this contribution, we focus on weighted null-space fitting
(WNSF), introduced in [39]. This method uses two of the
features of the methods above: i) an intermediate high-order
ARX model; ii) the high-order model is directly used for estimating the low-order model using ML-based model reduction.
However, instead of an explicit minimization of the modelreduction cost function—as in indirect PEM (directly via the
model parameters), ASYM (in the time domain), and [35]
(in the frequency domain)—the model reduction step consists
of a weighted least-squares problem. Asymptotic efficiency
requires that the weighting depend on the (to be estimated)
model parameters. To handle this, an additional least-squares
step is introduced. Consisting of three (weighted) least-squares
steps, WNSF has attractive computational properties in comparison with, for example, PEM, ASYM, and BJSM. More
steps may be added to this standard procedure, using an
iterative weighted least-squares algorithm, which may improve
the estimate for finite sample size.
Another interesting feature of WNSF is that, unlike many of
the methods above (including MORSM), the dynamic model
and the noise model are estimated jointly. If this is not done,
an algorithm cannot be asymptotically efficient for closedloop data [40]. Nevertheless, in some applications, the noise
model may be of no concern. WNSF can then be simplified
and a noise model not estimated, still maintaining asymptotic
efficiency for open-loop data. In closed loop, consistency
is still maintained because the high-order model captures
the noise spectrum consistently, while the resulting accuracy
corresponds to the covariance of PEM with an infinite-order
noise model [40]. Thus, besides the attractive numerical properties, WNSF has theoretical properties matched only by PEM.
However, WNSF has the additional benefit that an explicit
noise model is not required to obtain consistency with closedloop data.
In [39], some theoretical properties of WNSF are claimed
and supported by simulations, but with no formal proof. The
robust performance that the method has shown has provided
the motivation to extend the simulation study and deepen
the theoretical analysis. Take Fig. 1 as an example, showing
the FITs (see (V-C) for a definition of this quality measure)
for estimates obtained in closed loop from 100 Monte Carlo
runs with the following methods: PEM with default MATLAB
implementation (PEM-d), the subspace method SSARX [11],
WNSF, and PEM initialized at the true parameters as benchmark (PEM-t). Here, the default MATLAB initialization for
PEM is often not accurate enough, and the non-convex cost
function of PEM converges to non-global minima, while
the low FIT of SSARX indicates that this method is not
a viable alternative to deal with the non-convexity of PEM
for the situation at hand. On the other hand, WNSF has a

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

3

B. Definitions and Assumptions

FIT

100

Assumption 1 (Model and true system). The model has input
{ut }, output {yt } and is subject to the noise {et }, all realvalued, related by

50

yt = G(q, θ)ut + H(q, θ)et .
0

↓6

PEM-d

SSARX

WNSF

PEM-t

The transfer functions G(q, θ) and H(q, θ) are rational functions in q −1 , according to
L(q, θ)
l1 q −1 + · · · + lml q −ml
:=
,
F (q, θ)
1 + f1 q −1 + · · · + fmf q −mf
1 + c1 q −1 + · · · + cmc q −mc
C(q, θ)
:=
,
H(q, θ) :=
D(q, θ)
1 + d1 q −1 + · · · + dmd q −md

Fig. 1. FITs from 100 Monte Carlo runs with a highly resonant system.

performance close to PEM initialized at the true parameters,
suggesting that the weighted least-squares procedure applied
to a non-parametric estimate may be more robust than an
explicit minimization of the PEM cost function in regards to
convergence issues.
In this paper, we provide a theoretical and experimental
analysis of WNSF applied to stable single-input single-output
(SISO) Box-Jenkins (BJ) systems, which may operate in
closed loop. Our main contributions are to establish conditions
for consistency and asymptotic efficiency. A major effort of
the analysis is to keep track of the model errors induced
by using an ARX model on data generated by a system of
BJ type. It is a delicate matter to determine how the ARXmodel order should depend on the sample size such that it is
ensured that these errors vanish as the sample size grows: to
this end, the results in [41] have been instrumental. We finally
conduct a finite sample simulation study where WNSF shows
competitive performance with state-of-the-art methods.
The paper is organized as follows. In Section II, we introduce definitions, assumptions, and background. In Section III,
we review the WNSF algorithm. In Section IV, we provide the
theoretical analysis; in Section V, the experimental analysis.
II. P RELIMINARIES
A. Notation
•
•

•
•
•
•

•
•
•
•

(1)

G(q, θ) :=

where θ is the parameter vector to be estimated, given by

⊤
θ = f ⊤ l⊤ c⊤ d⊤ ∈ Rmf +ml +mc +md ,
(2)
with


f1


f =  ...  ,
f mf



l1
 
l =  ...  ,
l ml



c1


c =  ...  ,
cmc



d1


d =  ...  .
dmd


If the noise model is not of interest, we consider that we want
to obtain an estimate G(q, θ̄), where θ̄ = [f ⊤ l⊤ ]⊤ .
The true system is described by (1) when θ = θo . The
transfer functions Go := G(q, θo ) and Ho := H(q, θo ) are
assumed to be stable, and Ho inversely stable. The polynomials
Lo := L(q, θo ) and Fo := F (q, θo ), as well as Co := C(q, θo )
and Do := D(q, θo ), do not share common factors.
Because we allow for data to be collected in closed loop,
the input {ut } is allowed to have a stochastic part. Then, let
Ft−1 be the σ-algebra generated by {es , us , s ≤ t − 1}. For
the noise, the following assumption applies.
Assumption 2 (Noise). The noise sequence {et } is a stochastic process that satisfies

P
E[et |Ft−1 ] = 0, E[e2t |Ft−1 ] = σo2 , E[|et |10 ] ≤ C, ∀t.
||x||p = ( nk=1 |xk |p )1/p , with xk the k th entry of the
n × 1 vector x, and p ∈ N (for simplicity ||x|| := ||x||2 ).
Before stating the assumption on the input sequence, we
||A||p = supx6=0 ||Ax||p /||x||p , with A a matrix, x introduce the following definitions, used in [41].
a vector of appropriate dimensions, and p ∈ N (for
Definition 1 (fN -quasi-stationarity). Let fN be a decreasing
simplicity ||A|| := ||A||2 ); also, ||A||∞ = ||A⊤ ||1 .
C and N̄ denote any constant, which need not be the same sequence of positive scalars, with fN → 0 as N → ∞, and
 1 PN
in different expressions, and may be random variables.
⊤
,
0 ≤ τ < N,
 N t=τ +1 vt vt−τ
PN +τ
Γn (q) = [q −1 · · · q −n ]⊤ , where q −1 is the backward
N
1
⊤
Rvv (τ ) =
v
v
,
−N
< τ ≤ 0,
time-shift operator.
 N t=1 t t−τ
0,
otherwise.
∗
A is the complex conjugate transpose of the matrix A.
Tn,m (X(q)) is the Toeplitz matrix of size n × m (m ≤ The vector sequence {vt } is fN -quasi-stationary if
n) with first column [x0 · · ·P xn−1 ]⊤ and first row
i) There exists Rvv (τ ) such that
∞
[x0 01×m−1 ], where X(q) = k=0 xk q −k . The dimenN
sup|τ |≤N Rvv
(τ ) − Rvv (τ ) ≤ C1 fN ,
P
sion n may be infinity, denoted T∞,m (X(q)).
N
2
1
ii) N t=−N kvt k ≤ C2
Ex denotes expectation
of the random vector x.
P
N
for all N large enough, where C1 and C2 are finite constants.
Ēxt := lim N1 t=1 Ext .
N →∞
xN = O(fN ): the function xN tends to zero at a rate not
This definition allows us to work with some stochastic
slower than fN , as N → ∞, w.p.1.
signals that have deterministic components, as in [1]. In
xN ∼ AsN (a, P ): the random variable xN is normally addition to the standard definition of quasi-stationarity, a rate
distributed with mean a and covariance P as N → ∞.
of convergence for the sample covariances is defined.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

4

Definition 2P
(fN -stability). A filter G(q) =
∞
fN -stable if k=0 |gk |/fk < ∞.

P∞

gk q −k is

In open loop, the asymptotic covariance of the dynamicmodel parameters is the top-left block of (II-C) even if the
noise-model orders mc and md are larger than the true ones;
Definition 3 (Power spectral density). The power spectral
if smaller, the dynamic-model estimates are consistent but
density ofPan fN -quasi-stationary sequence {vt } is given by
not asymptotically efficient. In closed loop, the covariance of
∞
Φv (z) = τ =−∞ Rvv (τ )z −τ , if the sum exists for |z| = 1.
the dynamic-model estimates only corresponds to the top-left
For the input, the following assumption applies.
block of (II-C) if the noise-model orders are the true ones;
if smaller, the dynamic-model estimates are biased; if larger,
Assumption 3 (Input). The input sequence {ut } is defined by
they are consistent and the asymptotic covariance matrix can
ut = −K(q)yt + rt under the following conditions.
−1
be bounded by σo2 MCL
, where [40]
i) The sequence {rt } p
is independent of {et }, fN -quasiZ π
1
stationary with fN = log N/N , and uniformly bounded.
Ω̄(eiω )Φru (eiω )Ω̄∗ (eiω )dω
(7)
MCL =
−
1
2π
ii) With Φr (z) = Fr (z)Fr (z ) the spectral factorization of
−π
{rt } and Fr (z) causal, Fr (q) is BIBO stable.
with Φru the spectrum of the input due to the reference. This
√
iii) The closed loop system is fN -stable with fN = 1/ N .
corresponds to the case with infinite noise-model order.
iv) The transfer function K(z) is bounded on the unit circle.
The main drawback with PEM is that minimizing (II-C)
v) The spectral density of {[rt et ]⊤ } is coercive (i.e., is in general a non-convex optimization problem. Therefore,
PEM
bounded from below by the matrix δI, for some δ > 0). the global minimizer θ̂N
is not guaranteed to be found. An
exception
is
the
ARX
model.
Operation in open loop is obtained by taking K(q) = 0.
k=0

The choice of fN in iii) guarantees that the impulse responses
of the closed-loop system have a minimum rate of decay,
necessary to derive the results in [41]. This minimum decay
rate is trivially satisfied here, as the system is stable and finite
dimensional, and hence has exponentially decaying impulse
responses.

D. High-Order ARX Modeling
The true system can alternatively be written as
Ao (q)yt = Bo (q)ut + et ,

(8)

where the transfer functions
∞

X
1
=: 1+
aok q −k ,
Ao (q) :=
Ho (q)

C. The Prediction Error Method
The prediction error method minimizes a cost function of
prediction errors, which, for the model structure (1), are


L(q, θ)
D(q, θ)
yt −
ut .
εt (θ) =
C(q, θ)
F (q, θ)
Using a quadratic cost function, the PEM estimate of θ is
obtained by minimizing
PN 2
(3)
J(θ) = N1
t=1 εt (θ),

where N is the sample size. Assuming that θ belongs to
an appropriate domain [1, Def. 4.3], when the data set is
informative [1, Def. 8.1] and under appropriate technical
PEM
conditions [1, Chap. 8], the global minimizer θ̂N
of (II-C)
is asymptotically distributed as [1, Theorem 9.1]
√

N (θ̂PEM − θ ) ∼ AsN 0, σ 2 M −1 ,
(4)
o

N

o

CR

where

MCR

1
=
2π

Z

π

Ω(eiω )Φz (eiω )Ω∗ (eiω )dω,

(5)

k=1

A(q, η n )yt = B(q, η n )ut + et ,
where

η n = a1

A(q, η n ) = 1 +

· · · an

n
X

· · · bn

b1

ak q − k ,

B(q, η n ) =

k=1

(10)
⊤

,

n
X

bk q − k ,

k=1

can approximate (II-D) arbitrarily well if the model order n
is chosen arbitrarily large.
Because the prediction errors for the ARX model (II-D),
εt (η n ) = A(q, η n )yt − B(q, η n )ut , are linear in the model
parameters η n , the corresponding PEM cost function (II-C)
can be minimized with least squares. This is done as follows.
First, re-write (II-D) in regression form as
yt = (ϕnt )⊤ η n + et ,

iω

⊤

(9)

bok q −k

are stable (Assumption 1). Therefore, the ARX model

−π

with (for notational simplicity, we omit the argument e )

 Go
0
− Ho Fo Γmf
1


0
Ho Fo Γml

(6)
Ω=
1


0
C o Γmc
0
− D1o Γmd

k=1
∞
X

Go (q)
=:
Bo (q) :=
Ho (q)

(11)

where

ϕnt = −yt−1

· · · −yt−n

ut−1

· · · ut−n

⊤

.

Then, the least-squares estimate of η n is obtained by
η̂ n,ls = [Rn ]−1 rn ,
N

and Φz the spectrum of [ut et ] . When the error sequence is where
Gaussian, PEM with a quadratic cost function is asymptotiN
1 X n n ⊤
cally efficient, with (II-C) corresponding to the CR bound [1,
n
ϕ (ϕ ) ,
RN
=
N t=n+1 t t
Chap. 9].

N

N

n
rN
=

N
1 X n
ϕ yt .
N t=n+1 t

(12)

(13)

(14)

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

As the sample size increases, we have [41]


n
RN
→ R̄n := Ē ϕnt (ϕnt )⊤ , as N → ∞, w.p.1,

n
as N → ∞, w.p.1.
rN
→ r̄n := Ē [ϕnt yt ] ,

5

Then, (III) can be expanded as
(15)

∞


X
aok q −k
(1 + co1 q −1 + · · · + comc q −mc ) 1 +
k=1

Consequently,
n,ls
η̂N

 −1
→ η̄ := R̄n r̄n , as N → ∞, w.p.1.
n

For future reference, we define

⊤
ηon := ao1 · · · aon bo1 · · · bon ,

⊤
ηo := ao1 ao2 · · · bo1 bo2 · · · .

(19a)

− (1 + do1 q −1 + · · · + domd q −md ) = 0,

(16)

o
(1 + f1o q −1 + · · · + fm
q −mf )
f

−

(l1o q −1

+ ···+

∞
X

bok q −k

k=1
(19b)
∞


X
1+
aok q −k = 0.

o
lm
q −ml )
l

k=1

(17)

To express θo in terms of ηo , we can
so in vector
P∞form. BePdo
∞
cause a power-series product (e.g., k=0 αk q −k k=0 βk q −k )
The attractiveness of ARX modeling is the simplicity of es- can be written as the (Toeplitz-)matrix-vector product

  
 
timation while approximating more general classes of systems
α0
β0
β0
α0
0
0
 β1   β1 β0
 α1 
with arbitrary accuracy. However, as the order n typically has α1 α0

  
 
 β2  = β2 β1 β0
 α2  ,
to be taken large, the estimated ARX model will have high α2 α1 α0

  
 
variance. Nevertheless, this estimate can in principle be used
..
..
.. . . . . . .
.. . . . . . .
.
.
.
.
.
.
.
.
.
.
as a means to obtain an asymptotically estimate of the Box(20)
Jenkins (BJ) model (1) when the measurement noise is Gauswe may write (1) as (keeping the first n equations)
sian. To understand intuitively why this is the case, we observe
that if we neglect the truncation error from approximating
ηon − Qn (ηon )θo = 0,
(21)
n,ls
n
the true system (II-D) by the model (II-D), η̂N
and RN
with θo defined by (1) evaluated at the true parameters and
constitute a sufficient statistic for our problem. Therefore, they


0
0
−Qcn (η n ) Qdn
can replace the data without loss of information. If ML is used
n
Qn (η ) =
, (22)
−Qfn (η n ) Qln (η n )
0
0
for the subsequent estimation, we need to solve a non-convex
optimization problem [35]. An accurate estimate to initialize where, when evaluated at the true parameters η n ,
o
the optimization procedure is then crucial; a standard result
n
l
c
n
Qn (ηo ) = Tn,mc (A(q, ηo )), Qn (ηo ) = Tn,ml (A(q, ηo )),
is that, if initialized with a strongly consistent estimate, one


Gauss-Newton iteration provides an asymptotically efficient
Imd ,md
d
f
n
Qn (ηo ) = Tn,mf (B(q, ηo )), Qn =
.
estimate (e.g., [42, Chap. 23]).
0n−md ,md
III. W EIGHTED N ULL -S PACE F ITTING M ETHOD
The idea of weighted null-space fitting [39] is to avoid
the burden of a non-convex optimization by using weighted
least squares, but maintaining the properties of maximum
likelihood. The method consists of three steps. In the first
step, a high-order ARX model is estimated with least squares.
In the second step, the parametric model is estimated from
the high-order ARX model with least squares, providing a
consistent estimate. In the third step, the parametric model is
re-estimated with weighted least squares. Because the optimal
weighting depends on the true parameters, we replace these by
the consistent estimate obtained in the previous step. Similarly
to maximum likelihood with an optimization algorithm initialized at a consistent estimate, this provides an asymptotically
efficient estimate. We now proceed to detail each of these
steps.
n,ls
The first step consists in estimating η̂N
from (II-D). As
n,ls
n
discussed before, η̂N and RN are almost a sufficient statistic
for our problem, if the ARX-model truncation error is small
enough (later, this will be treated formally). Then, we will use
n,ls
n
η̂N
and RN
instead of data to estimate the model of interest.
The second step implements this as follows. Re-write (II-D)
as
Co (q)Ao (q) − Do (q) = 0,
(18)
Fo (q)Bo (q) − Lo (q)Ao (q) = 0.

n,ls
Motivated by (III), we replace ηon by its estimate η̂N
, obtaining an over-determined system of equations, which may be
solved for θ using, for example, least squares:
−1

n,ls n,ls
n,ls
n,ls
LS
Q⊤
(23)
)
)Q
(η̂
θ̂N
= Q⊤
(η̂
n
n (η̂N )η̂N .
n N
N

n,ls
In (III), invertibility follows from convergence of η̂N
to ηon ,
which is of larger dimension than θ, and the block-Toeplitz
structure of Qn (η n ) (this is treated formally in Lemma 1).
With (III), we have not accounted for the residuals in (III)
n,ls
when η̂N
replaces ηon . The third step remedies this by reestimating θ in a statistically sound way.
For some η n , and using the same logic as (III), we can
write (III) as

η n − Qn (η n )θo = Tn (θo )(η n − ηon ) =: δn (η n , θo ).
where

Tnc (θ)
Tn (θ) =
−Tnl (θ)



0
,
Tnf (θ)

(24)
(25)

with Tnc (θo ) = Tn,n (C(q, θo )), Tnl (θo ) = Tn,n (L(q, θo )), and
Tnf (θo ) = Tn,n (F (q, θo )). The objective is then to estimate
n,ls
θ that minimizes the residuals δn (η̂N
, θ). If we neglect
the bias error from truncation of the ARX model, which
should be close to zero for sufficiently large n, we have that,
approximately,
√
 −1 

n,ls
N η̂N
− ηon ∼ AsN 0, σo2 R̄n
.
(26)

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

6

Then, using (III) and (III), we may write that, approximately,


n,ls
δn (η̂N
) ∼ AsN 0, Tn (θo )σo2 [R̄n ]−1 Tn⊤ (θo ) .
(27)
Because the residuals we want to minimize, given by
n,ls
n,ls
n,ls
δn (η̂N
, θ) = η̂N
−Qn (η̂N
)θ, are asymptotically distributed
by (III), the estimate of θ with minimum variance is given by
the weighted least-squares estimate

−1
n,ls
n,ls
n,ls
n,ls
WLSo
Q⊤
θ̂N
= Q⊤
n (η̂N )W̄n (θo )η̂N ,
n (η̂N )W̄n (θo )Qn (η̂N )

where the weighting matrix
−1

W̄n (θo ) = Tn (θo )σo2 [R̄n ]−1 Tn⊤ (θo )

(28)

is the inverse of the covariance of the residuals [43]. Because
LS
n
θo and R̄n are not available, we replace them by θ̂N
and RN
,
2
respectively (σo can be disregarded, because the weighting can
be scaled arbitrarily without influencing the solution). Thus,
the third step consists in re-estimating θ by

−1
n,ls
n,ls
n,ls
WLS
LS
LS n,ls
θ̂N
= Q⊤
(η̂
)W
(
θ̂
)Q
(η̂
)
Q⊤
n
n
n
N
n (η̂N )Wn (θ̂N )η̂N ,
N
N
(29)
where (we take the inverses of the matrices individually)
W (θ̂LS ) = T −⊤ (θ̂LS )Rn T −1 (θ̂LS ),
(30)
n

N

n

N

N

n

N

LS
Tn (θ̂N
)

with
obtained using (III). Invertibility in (III) follows
from (besides what was mentioned for Step 2) invertibility
LS
of Wn (θ̂N
), which in turn follows from the lower-Toeplitz
LS
n
structure of Tn (θ), convergence of θ̂N
to θo and RN
to
n
R̄ (this is treated formally in Lemmas 2 and 3). Because
LS
θ̂N
is a consistent estimate of θo with an error decaying
LS
sufficiently fast, using θ̂N
in the weighting should not change
WLS
the asymptotic properties of θ̂N
. Analogously to taking
one Gauss-Newton iteration for a maximum likelihood cost
WLS
function initialized at a strongly consistent estimate, θ̂N
is
an asymptotically efficient estimate, as will be proven in the
next section.
In summary, WNSF consists of the following three steps:
1) estimate a high-order ARX model with least
squares (II-D);
2) reduce the high-order ARX model to the model of interest
with least squares (III);
3) re-estimate the model of interest by weighted least
squares (III) using the weighting (III).
Two notes can be made about this procedure. First, the
objective of the second step is to obtain a consistent estimate
to construct the weighting; hence, the choice of least squares
is arbitrary, and weighted least squares with any invertible
n
weighting (e.g., Wn = RN
) can be used. Second, although
WLS
θ̂N is asymptotically efficient, it is possible to continue
iterating, which may improve the estimate for finite sample
size.
Other Settings
Despite having been presented for a fully parametrized
SISO BJ model, we point out that the method is flexible
in parametrization. For example, it is possible to fix some

parameters in θ if they are known, or to impose linear relations
between parameters. Hence, other common model structures
(e.g., OE, ARMA, ARMAX) may also be used, as well as
multi-input multi-output (MIMO) versions of such structures.
The requirement is that a relation between the high- and loworder parameters can be written in the form (III).
Moreover, a parametric noise model does not need to be
estimated. In this case, disregard (1) and consider only (1). The
subsequent steps can then be derived similarly. This approach
is presented in detail and analyzed in [44]. In open loop,
it provides asymptotically efficient estimates of the dynamic
model; in closed loop, the estimates are consistent and with
asymptotic covariance corresponding to (II-C).
IV. A SYMPTOTIC P ROPERTIES
We now turn to the asymptotic analysis of WNSF. Here, we
make a distinction between the main algorithm presented here
and the aforementioned case without a low-order noise model
estimate. Although apparently simpler because of the smaller
dimension of the problem, the case without a noise-model
estimate requires additional care in the analysis. The reason is
that the corresponding Tn (θ) in that case will not be square.
Then, inverting the weighting as in (III) (a relation that will be
used for the analysis in this paper) will not be valid, requiring
another approach. Including in this paper under-parametrized
noise models is then not possible for space concerns. Thus,
the asymptotic analysis in this paper considers the dynamic
and noise models correctly parametrized, in which case the
algorithm is consistent and asymptotically efficient. The case
with an under-parametrized noise model (in particular, the
limit case where (1) is neglected and no noise-model is
estimated) is considered in [44].
Because the ARX model (II-D) is a truncation of the true
system (II-D), its estimate (and the respective covariance)
will not be a sufficient statistic for finite order, and some
information will be lost in this step. Then, we need to make
sure that, as N grows, the truncation error will be sufficiently
small so that, asymptotically, no information is lost. To keep
track of the truncation error in the analysis (see appendices),
we let the model order n depend on the sample size N —
denoting n = n(N )—according to the following assumption.
Assumption 4 (ARX-model order). It holds that
D1. n(N ) → ∞, as N → ∞;
D2. n4+δ (N )/N → 0, for some δ > 0, as N → ∞.
Condition D1 implies that, as the sample size N tends to
infinity, so does the model order n. Condition D2 establishes a
maximum rate at which the model order n is allowed to grow,
as we cannot use too high order compared with the number
of observations. A consequence of Condition D2 is that [41]
n2 (N ) log(N )/N → 0, as N → ∞,
(31)
P∞
Moreover, defining d(N ) := k=n(N )+1 |aok | + |bok | , we have
√
N d(N ) → 0, as N → ∞,
(32)
as consequence of stability and rational description of the true
system in Assumption 1. Although (IV) and (IV) follow from

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

7

other assumptions, they are stated explicitly as they will be
required to show our theoretical results.
To facilitate the statistical analysis, the results in this section
consider, instead of (II-D), a regularized estimate
n,reg
n
n
n
η̂N
:= η̂N
= [Rreg
(N )]−1 rN
,

(33)

where
n
Rreg
(N )

=



n
RN
n
RN
+ 2δ I2n

n −1
if ||[RN
] || < 2/δ,
otherwise,

for some small δ > 0. Asymptotically, the first and second
n,ls
n
order properties of η̂N
and η̂N
are identical [41].
When we let n = n(N ) according to Assumption 4, we
n(N )
n(N )
,
use η̂N := η̂N . We will also denote η̄ n(N ) and ηo
defined in (II-D) and (II-D), respectively. Concerning the
matrices (II-D), (III), (III), (III), and (III), for notational
simplicity we maintain the subscript n even if n = n(N ).
Some of the technical assumptions used in this paper differ
from those used for the asymptotic analysis of PEM [1]. For
example, the bound in Assumption 2 is stronger than what is
required for PEM. On the other hand, for PEM the parameter
vector θ is required to belong to a compact set, which is
not imposed here. However, such differences in technical
assumptions have little relevance in practice.
LS
We have the following result for consistency θ̂N
.
LS
Theorem 1. Let Assumptions 1, 2, 3, and 4 hold, and and θ̂N
be defined by (III). Then,
LS
θ̂N

In this section, we perform simulation studies and discuss
practical issues. First, we illustrate the asymptotic properties
of the method. Second, we consider how to choose the order
of the non-parametric model. Third, we exemplify with two
difficult scenarios for PEM how WNSF can be advantageous in
terms of robustness against convergence to non-global minima
and convergence speed. Fourth, we perform a simulation with
random systems to test the robustness of the method compared
with other state-of-the-art methods.
Although WNSF and the approach in [38] are different
algorithms, they share the similarities of using high-order
models and iterative least squares. However, [38] is only
applicable in open loop. Here, to differentiate WNSF as a
more general approach that is applicable in open or closed
loop without changing the algorithm, we focus on the typically
more challenging closed-loop setting, for which many standard
methods are not consistent.
A. Illustration of Asymptotic Properties
The first simulation has the purpose of illustrating that the
method is asymptotically efficient. Here, we consider only
the case where we estimate a correct noise model (the case
where a low-order noise model is not estimated is illustrated in
[44]). We perform open- and closed-loop simulations, where
the closed-loop data are generated by
K(q)Ho (q)
1
rt −
et ,
1 + K(q)Go (q)
1 + K(q)Go (q)
Go (q)
Ho (q)
yt =
rt +
et ,
1 + K(q)Go (q)
1 + K(q)Go (q)

ut =

→ θo , as N → ∞, w.p.1.

Moreover, we have that

r

log N
LS
1 + d(N ) .
n(N )
||θ̂N − θo || = O
N

V. S IMULATION S TUDIES

(34)

Proof. See Appendix B.

WLS
We have the following result for consistency of θ̂N
.
WLS
Theorem 2. Let Assumptions 1, 2, 3, and 4 hold, and θ̂N
be defined by (III). Then,
WLS
θ̂N
→ θo , as N → ∞, w.p.1.

Proof. See Appendix C.
We have the following result for asymptotic distribution and
WLS
covariance of θ̂N
.
WLS
Theorem 3. Let Assumptions 1, 2, 3, and 4 hold, and θ̂N
be defined by (III). Then,
√
−1 ),
WLS
N (θ̂N
− θo ) ∼ AsN (0, σo2 MCR

where MCR is given by (II-C).
Proof. See Appendix D.
Theorem 3 implies, comparing with (II-C), that WNSF has
the same asymptotic properties as PEM. For Gaussian noise,
this corresponds to an asymptotically efficient estimate.

(35)

and the open-loop data by
ut =

1
rt ,
1 + K(q)Go (q)

yt = Go (q)ut + Ho (q)et ,

where {rt } and {et } are independent Gaussian white sequences with unit variance, K(q) = 1, and
Go (q) =

q −1 + 0.1q −2
,
1 − 0.5q −1 + 0.75q −2

Ho (q) =

1 + 0.7q −1
.
1 − 0.9q −1

We perform 1000 Monte Carlo runs, with sample sizes
N ∈ {300, 600, 1000, 3000, 6000, 10000}. We apply WNSF
with an ARX model of order 50 with the open- and closedloop data. Performance is evaluated by the mean-squared
error of the estimated parameter vector of the dynamic model,
WLS
− θ̄o ||2 , where θ̄ contains only the elements of
MSE = ||θ̄ˆN
θ contributing to G(q, θ). As this simulation has the purpose
of illustrating asymptotic properties, initial conditions are zero
and assumed known—that is, the sums in (II-D) start at t = 1
instead of t = n + 1.
The results are presented in Fig. 2, with the average MSE
over 1000 Monte Carlo runs plotted as function of the sample
size (closed loop in solid line, open loop in dash-dotted line),
where we also plot the corresponding CR bounds (closed loop
in dashed line, open loop in dotted line). The respective CR
bounds are attained as the sample size increases.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

8

TABLE I
C OMPARISON WITH PEM: AVERAGE FIT S WITH DIFFERENT METHODS
(M ETH ) AND INITIALIZATIONS (I NIT ).

MSE

10−2

Init

10−3

10−4
103

N

104

Fig. 2. Illustration of asymptotic properties: CR bounds in closed loop
(dashed) and open loop (dotted), and average MSE for the dynamic-model
parameter estimates as function of sample size obtained with WNSF in closed
loop (solid) and open loop (dash-dotted).

B. Practical Issues

Meth
WNSF
PEM GN
PEM LM

MtL

LS

true

98
74
98

98
87
85

–
98
98

with {et } and {rtw } Gaussian white noise sequences with variances 4 and 0.25, respectively. The sample size is N = 2000.
We estimate an OE model with the following algorithms:
• WNSF with a non-parametric model of order n = 250;
• PEM with default MATLAB initialization and GaussNewton (GN) algorithm;
• PEM with default MATLAB initialization (MtL) and
Levenberg-Marquardt (LM) algorithm;
• WNSF with a non-parametric model of order n = 250,
where the weighting matrix, instead of being initialized
LS
with θ̂N
(III), is initialized with the default MATLAB
initialization (MtL);
LS
• PEM initialized with θ̂N (LS) and the GN algorithm;
LS
• PEM initialized with θ̂N (LS) and the LM algorithm;
• PEM initialized at the true parameters (true).
All the methods use a maximum of 100 iterations, but stop
early upon convergence (default settings for PEM, 10−4 as
tolerance for the normalized relative change in the parameter
estimates) and initial conditions are zero.
Performance is evaluated by the FIT of the impulse response
WLS
of the estimated OE model G(q, θ̂N
), given in percent by


kgo − ĝk
,
(36)
FIT = 100 1 −
kgo − mean(go )k

In the previous simulation, an ARX model of order 50
was estimated in the first step. Although the order of this
model should, in theory, tend to infinity at some maximum
rate to attain efficiency (Assumption 4), a fixed order was
sufficient to illustrate the asymptotic properties of WNSF in
this particular scenario. This suggests that when the number of
data samples increases, a non-parametric model of fixed order
with sufficiently low bias error may be enough for practical
purposes. However, for fixed sample size, the question remains
on how to choose the most appropriate non-parametric model
order: a too small n will introduce bias, and a too large n
will introduce unnecessary noise in the non-parametric model
estimate, which may affect the accuracy of the parametric
model estimate. Some previous knowledge about the speed
of the system may help in choosing this order, but the
most appropriate value may also depend on sample size and
signal-to-noise ratio. In this paper, we use the PEM cost
WLS
function (II-C) as criterion to choose n: we compute θ̂N
for several n, and choose the estimate that minimizes (II-C). where go is a vector with the impulse response parameters
WLS
need not be used as final estimate, as, for finite of Go (q), and similarly for ĝ but for the estimated model.
Also, θ̂N
sample size, performance may improve by iterating. However, In (V-C), sufficiently long impulse responses are taken to make
because WNSF does not minimize the cost function (II-C) sure that the truncation of their tails does not affect the FIT.
The average FITs for 100 Monte Carlo runs are shown
explicitly, it is not guaranteed that subsequent iterations correin
Table I. For PEM, the results depend on the optimization
spond to a lower cost-function value than previous ones. Here,
method
and the initialization point: as consequence of the nonwe will also use the cost function (II-C) as criterion to choose
convexity
of PEM, the algorithms do not always converge to
the best model among the iterations performed.
the global optimum. For PEM implementations, the average
FIT is the same as for PEM started at the true parameters only
C. Comparison with PEM
with default MATLAB initialization and LM algorithm. For
One of the main limitations of PEM is the non-convex WNSF, the average FIT is the same as for PEM started at the
cost function, which may make the method sensitive to the true parameters independently of the initialization point used
initialization point. Here, we provide examples illustrating how in the weighting matrix, suggesting robustness to different
WNSF may be a more robust method than PEM regarding initial weighting matrices.
initialization: in cases where the PEM cost function is highly
In this simulation, PEM was most robust with the LM
non-convex, WNSF may require less iterations and be more algorithm and the default MATLAB initialization, having on
robust against convergence to non-global minima.
average the same accuracy as WNSF. Then, it is appropriate to
We consider a system where Ho (q) = 1, K(q) = 0.3, and
compare the performance of these methods by iteration when
WNSF
also is initialized with the same parameter values. In
−1
−2
1.0q − 1.2q
,
Go (q) =
Fig.
3,
we
plot the average FITs for these methods as function
1 − 2.5q −1 + 2.4q −2 − 0.88q −3
of the maximum number of iterations. Here, WNSF reaches
with data generated according to (V-A), where
an average FIT of 98 after two iterations, while PEM with
−1
LM
takes 20 iterations to reach the same value. This suggests
1 + 0.7q
rt =
rtw ,
that,
even if WNSF and some PEM implementation start and
−1
1 − 0.9q

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

9

For the simulation, we use 100 systems with structure

100

average FIT

80

Go (q) =

60
40
20

WNSF
PEM LM

0
−20
0

5

10

15

max iterations

20

25

Fig. 3. Comparison with PEM: average FIT from 100 Monte Carlo runs
function of the maximum number of iterations.

converge to the same value, WNSF may do it faster than
standard optimization methods for PEM.
The robustness of WNSF against convergence to non-global
minima compared with different instances of PEM can be even
more evident than in Table I, as WNSF seems to be appropriate
for modeling systems with many resonant peaks, for which the
PEM cost function can be highly non-linear. Take the example
in Fig. 1, based on 100 Monte Carlo runs for a system with
Lo (q) = q −1 − 3.4q −2 + 4.8q −3 − 3.3q −4 + 0.96q −5 ,

As we have observed, PEM may have difficulties with slow
resonant systems: therefore, it is for this class of systems that
WNSF may be most beneficial. With this purpose, we generate
the polynomial coefficients in the following way. The poles are
located in an annulus with the radius uniformly distributed
between 0.88 and 0.98, and the phase uniformly distributed
between 0 and 90◦ (and respective complex conjugates). One
pair of zeros is generated in the same way, and a third real
zero is uniformly distributed between −1.2 and 1.2 (this
allows for non-minimum-phase systems). The noise models
have structure
1 + co1 q −1 + co2 q −2
,
Ho (q) =
1 + do1 q −1 + do2 q −2
with the poles and zeros having uniformly distributed radius
between 0 and 0.95, and uniformly distributed phase between
0 and 180◦ (and respective complex conjugates).
The data are generated in closed loop by

Fo (q) = 1 − 5.4q −1 + 13.5q −2 − 20.1q −3 + 19.5q −4
− 12.1q −5 + 4.5q −6 ,

and data generated according to (V-A) with K(q) = −0.05,
0.05
rt =
rw ,
1 − 0.99q −1 t
{rtw }

where
and {et } are Gaussian white sequences with
unit variance. Here, initial conditions are not assumed zero:
PEM estimates initial conditions by backcasting and WNSF
uses the approach in [45]. In this scenario, PEM with the
LM algorithm and default initialization fails in most runs
to find the global optimum. Subspace methods, often used
to avoid the non-convexity of PEM, may not help in this
scenario: SSARX [11], a subspace method that is consistent
in closed loop, provides an average FIT around 20% (default
MATLAB implementation). Here, WNSF with n between 100
and 600 spaced with intervals of 50 performs similarly to
PEM initialized at the true parameters, accurately capturing
the resonance peaks of the system.

l1o q −1 + · · · + l4o q −4
.
o q −6
1 + f6o q −6 + · · · + fm

K(q)
K(q)Ho (q)
rt −
et ,
1 + K(q)Go (q)
1 + K(q)Go (q)
Ho (q)
K(q)Go (q)
rt +
et ,
yt =
1 + K(q)Go (q)
1 + K(q)Go (q)

ut =

where
rt =

1 − 1.273q −1 + 0.81q −2 w
r
1 − 1.559q −1 + 0.81q −2 t

with {rtw } a Gaussian white-noise sequence with unit variance, {et } a Gaussian white-noise sequence with the variance
chosen such that the signal-to-noise ratio (SNR) is
PN h K(q)Go (q) i2
t=1 1+K(q)Go (q) rt
SNR =
= 2,
PN
2
t=1 [Ho (q)et ]

and the controller K(q) is obtained using a Youlaparametrization to have an integrator and a closed-loop transfer
function that has the same poles as the open loop except that
the radius of the slowest open-loop pole pair is reduced by
80%. The sample size is N = 2000 and we perform 100 Monte
Carlo runs (one for each system; different noise realizations).
D. Random Systems
We compare the following methods:
In order to test the robustness of the method, we now
• PEM initialized at the true parameters (PEMt);
perform a simulation with random systems. Also, closed-loop
• PEM with default MATLAB initialization (PEMd);
data often introduces additional difficulties: for example, many
• SSARX with the default MATLAB options;
standard methods are not consistent. Thus, we perform a
• WNSF using the approach in Section V-B to choose n
simulation with these settings and compare the performance of
from the grid {50, 100, 150, 200, 250, 300}.
WNSF with other methods available in the Mathworks MAT• PEM initialized with WNSF (PEMw).
LAB System Identification Toolbox. For a fair comparison, we
only use methods that are consistent in closed loop and only All methods estimate a fully parametrized noise model. We
use input and output data. From the subspace class, we use use the MATLAB2016b System Identification Toolbox impleSSARX, as this method is competitive with other subspace mentation of SSARX and PEM. For PEM, the optimization
algorithms such as CVA [46,47] and N4SID [48], while it is algorithm is LM. For SSARX, the horizons are chosen auconsistent in closed loop [11]. IV methods are not considered, tomatically by MATLAB, based on the Akaike Information
as they in closed loop require the reference signal to construct Criterion. WNSF and PEM use a maximum of 100 iterations,
but stop earlier upon convergence (default settings for PEM,
the instruments.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

10

Future work includes also extensions to dynamic networks and
non-linear model structures.

100

FIT

75
50

A PPENDIX A
AUXILIARY R ESULTS

25
0
−25

In this appendix, we present some results that will be
applied in the remainder of the paper.

↓3

PEMd

SSARX WNSF

PEMw

PEMt

Fig. 4. Random systems: FITs from 100 Monte Carlo runs.

10−4 as tolerance for the normalized relative change in the
parameter estimates. PEM estimates initial conditions by backcasting and WNSF truncates them ([45] does not apply to BJ
models).
The FITs obtained in this simulation are presented in Fig. 4.
In this scenario, PEM with default MATLAB initialization
(PEMd) often fails to find a point close to the global optimum,
which can be concluded by comparison with PEM initialized at
the true parameters (PEMt). Also, SSARX is not an alternative
for achieving better performance. WNSF can be an appropriate
alternative, failing only once to provide an acceptable estimate,
and having otherwise a performance close to the practically
infeasible PEMt. The estimate obtained with WNSF may be
used to initialize PEM. This provides a small improvement
only, suggesting that the estimates obtained with WNSF are
already close to a (local) minimum of the PEM cost function.

Proposition 1. Let Assumptions 1, 2, and 3 hold. Also, let η̄ n
be defined by (II-D) and ηon by (II-D). Then,
P∞
kη̄ n − ηon k ≤ C k=n+1 |aok | + |bok | → 0, as n → ∞.
Proof. The result follows from [41, Lemma 5.1] and (IV).

Proposition 2. Let Assumptions 1, 2, 3, and 4 hold. Also, let
n(N )
η̂N := η̂N
be defined by (IV) and η̄ n(N ) by (II-D). Then,
r

n(N ) log N
n(N )
η̂N − η̄
=O
[1 + d(N )] ,
(37)
N
and η̂N − η̄ n(N ) → 0, as N → ∞, w.p.1.
Proof. For the first part, see [41, Theorem 5.1]. The second
part follows from (IV) and (IV).
Proposition 3. Let Assumptions 1, 2, 3, and 4 hold. Then,
!
r
2
n
(N
)
log
N
n(N )
.
+C
RN − R̄n(N ) = O 2n(N )
N
N
Proof. See [41, Lemma 4.1].

VI. C ONCLUSION
Methods for parameter estimation based on an intermediate
unstructured model have a long history in system identification (e.g., [11,35]–[37]). Here, we believe to have taken a
significant step further in this class of methods, with a method
that is flexible in parametrization and provides consistent and
asymptotically efficient estimates in open and closed loop
without using a non-convex optimization or iterations.
In this paper, we provided a theoretical and experimental
analysis of this method, named weighted null-space fitting
(WNSF). Theoretically, we showed that the method is consistent and asymptotically efficient for stable Box-Jenkins systems. Experimentally, we performed Monte Carlo simulations,
comparing PEM, subspace, and WNSF under settings where
PEM typically performs poorly. The simulations suggest that
WNSF is competitive with these methods, being a viable
alternative to PEM or to provide initialization points for PEM.
Although WNSF was here presented for SISO BJ models,
it was also pointed out that the flexibility in parametrization
allows for a wider range of structures to be used, as well as
for incorporating structural information (e.g., fixing specified
parameters). Moreover, based on the analysis in [44], WNSF
does not require a parametric noise model to achieve asymptotic efficiency in open loop and consistency in closed loop.
An extension that was not covered in this paper is the
MIMO case, where subspace or IV methods are typically
used [49], as PEM often has difficulty with estimation of such
systems. Based on the theoretical foundation provided in this
contribution, this important extension is already in preparation.

Proposition 4. Let Assumptions 1, 2, 3, and 4 hold. Also, let
Υn be an m × 2n deterministic matrix, with m fixed. Then,
√
N Υn (η̂N − η̄ n(N ) ) ∼ AsN (0, P ),
where P = σo2 lim Υn [R̄n ]−1 (Υn )⊤ , if the limit exists.
n→∞

Proof. See [41, Theorem 7.3].
Qp
(i)
Proposition 5. Consider the product i=1 XN , where p
(i)
is finite and XN are stochastic matrices of appropriate
dimensions (possibly a function of N ) such that
(i)

||XN − X̄ (i) || → 0, as N → ∞, w.p.1
and ||X̄ (i) || < Ci . Then, we have that
Qp
Qp
(i)
(i)
→ 0, as N → ∞, w.p.1. (38)
i=1 X̄
i=1 XN −

Proof. We show this by induction. First, let p = 2 and define
(i)
(i)
(i)
∆N := XN − X̄N . Then, we can write
(1)

(2)

(1)

(2)

(1)

(2)

XN XN − X̄ (1) X̄ (2) = ∆N X̄ (2) + X̄ (1) ∆N + ∆N ∆N ,
(39)
which yields, using the assumptions,
(1)

(2)

(1)

||XN XN − X̄ (1) X̄ (2) || ≤ ||∆N || ||X̄ (2) ||
(2)
(1)
(2)
+||X̄ (1) || ||∆N ||+||∆N || ||∆N || → 0 as N → ∞, w.p.1.
(40)
Second, we consider an arbitrary p, and assume that
Qp−1 (i) Qp−1 (i)
→ 0, as N → ∞, w.p.1. (41)
i=1 X̄
i=1 XN −

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

11

Then, using a similar procedure as (A), we have
Qp−1
Qp
Qp
(p)
(i)
(i)
≤ ||∆N || || i=1 X̄ (i) ||
i=1 X̄
i=1 XN − Q
Qp−1 (i)
(i)
(p)
+ ||X̄ (p) || || p−1
i=1 ∆N || + ||∆N || ||
i=1 ∆N ||,
(42)
which, in turn, is bounded by
Q
(p) Q
(i)
(i)
(p)
||∆N || p−1
|| p−1
i=1 ||X̄ || + ||X̄
i=1 ||∆N ||
(i)
(p) Qp−1
+ ||∆N || i=1 ||∆N || → 0 as N → ∞, w.p.1,
where the convergence follows by assumption. Then, (5) is
verified when assuming (A), which considering also (A) and
an induction argument, concludes the proof.

• kTn (θo )k is bounded for all n
Consider the matrix Tn (θo ), givenP
by (III). First, we introduce
∞
the following result. Let X(q) = k=0 xk q −k and define


x0 0
0 ···


x1 x0 0 . . . 


T[X(q)] := 
(50)
.
x2 x1 x0 . . . 


.. . . . .
..
.
.
.
.
pP∞
2
If
k=0 |xk | < C1 , we have that [50]

kT[X(q)]k ≤ C.

(51)

A PPENDIX B
When X(q) can be written as a rational transfer function, (B)
C ONSISTENCY OF S TEP 2
follows from X(q) having all poles strictly inside the unit
The main purpose of this appendix is to prove Theorem 1. circle, as, in this case, the sum of squares of its impulse
However, before we do so, we introduce some results regard- response coefficients is bounded.
We observe that the blocks of Tn (θo ) satisfy that Tnf (θo ),
ing the norm of some vectors and matrices.
c
Tn (θo ), and Tnl (θo ) are sub-matrices of T[Fo (q)], T[Co (q)],
n(N )
•
η̂N − ηo
tends to zero, as N tends to infinity, w.p.1
and T[Lo (q)], respectively. Then, we have that
n(N )
Consider the estimated parameter vector η̂N := η̂N
(IV),
n(N )
kTn (θo )k ≤ Tnf (θo ) + kTnc (θo )k + Tnl (θo )
and the truncated true parameter vector ηo
(II-D). Using
≤ kT[Fo (q)]k + kT[Co (q)]k + kT[Lo (q)]k (52)
the triangular inequality, we have
η̂N − ηon(N ) ≤ η̂N − η̄ n(N ) + η̄ n(N ) − ηon(N ) , (43)

where η̄ n is defined by (II-D). Then, from Proposition 1,
the second term on the right side of (B) tends to zero as
n(N ) → ∞. From Proposition 2, the first term on the right
side of (B) tends to zero, as N → ∞, w.p.1. Thus,
η̂N − ηon(N ) → 0, as N → ∞, w.p.1.

n(N )
• Qn (η̂N )−Qn (ηo
) tends to zero, as N
n(N )
Consider Qn (ηo
), given by (III) evaluated
n(N )
true parameter vector ηo
, and the matrix

(44)

to infinity, w.p.1

at the truncated
Qn (η̂N ), given
by (III) evaluated at the estimated parameters η̂N . We have
Qn (η̂N ) − Qn (ηon(N ) ) ≤
+ Qln (η̂N ) − Qln (ηon(N ) )
≤ C η̂N − ηon(N ) .

Qcn (η̂N ) − Qcn (ηon(N ) )
+ Qfn (η̂N ) − Qfn (ηon(N ) )
(45)

Then, using (B), we conclude that
Qn (η̂N ) − Qn (ηon(N ) ) → 0, as N → ∞, w.p.1.

(46)

Qn (ηon ) is bounded for all n
We have that
•

kQn (ηon )k ≤ kQcn (ηon )k+ Qln (ηon ) + Qfn (ηon ) + Qdn
≤ C kηon k + 1 ≤ C kηo k + 1, ∀n

(47)

which is bounded, by stability of the true system.
• kQn (η̂N )k is bounded for large N , w.p.1
Using the triangular inequality, we have
kQn (η̂N )k ≤ Qn (η̂N )−Qn (ηon(N ) ) + Qn (ηon(N ) ) . (48)
Using now (B) and (B), the first term on the right side of (B)
can be made arbitrarily small as N increases, while the second
term is bounded for all n(N ). Then, there exists N̄ such that
kQn (η̂N )k ≤ C,

∀N > N̄ .

(49)

≤ C ∀n,

where the last inequality follows from (B) and from F (q),
C(q), and L(q) being finite order polynomials.
The following lemma is useful for invertibility of the leastsquares problem (III).
Lemma 1. Let Assumption 1 hold and
n
n
M (ηo ) := lim Q⊤
n (ηo )Qn (ηo ),
n→∞

(53)

where Qn (ηon ) is given by (III) evaluated at ηon , defined
in (II-D). Then, M (ηo ) is invertible.
Proof. Firts, we observe that the limit in (1) is well defined,
because the entries of M (ηon ) := Q⊤ (ηon )Q(ηon ) are either
zero or sums with form
Pn
Pn
Pn
o o
o o
o o
k=1 bk bk+p ,
k=1 ak bk+p ,
k=1 ak ak+p ,

for some finite integers p, and the coefficients aok and bok are
stable sequences. Thus, these sums converge as n → ∞. For
simplicity of notation, let Q∞ (ηo ) := limn→∞ Qn (ηon ); that is,
Q∞ (ηo ) is block Toeplitz according to (III), with each block
having an infinite number of rows and given by
Qc∞ (ηo ) = T∞,mc (A(q, ηo )), Ql∞ (ηo ) = T∞,ml (A(q, ηo )),


Imd ,md
f
d
Q∞ (ηo ) = T∞,mf (B(q, ηo )), Q∞ =
.
0∞,md

We can then write M (ηo ) = Q⊤
∞ (ηo )Q∞ (ηo ). From this
factorization, we observe that M (ηo ) is singular if and only if
Q∞ (ηo ) has a non-trivial right null-space. Moreover, the block
anti-diagonal structure of Q∞ (ηo ) implies that Q∞ has full
column rank if and only if both matrices [−Qf∞ (ηo ) Ql∞ (ηo )]
and [−Qc∞ (ηo ) Qd∞ (ηo )] have full column rank. We proceed
by contradiction. Suppose that
 

 α
f
l
−Q∞ (ηo ) Q∞ (ηo )
= −Qf∞ (ηo )α + Ql∞ (ηo )β = 0,
β
(54)

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

12

where α and β are some vectors α = [α0 . . . αmf −1 ]⊤ and
β = [β0 . . . βml−1 ]⊤ . Then, (B) implies
B(q, ηo )α(q) = A(q, ηo )β(q) ⇔ L(q, θo )α(q) = F (q, θo )β(q),
Pmf −1
Pml −1 (55)
−k
−k
where α(q) =
α
q
and
β(q)
=
.
k
k=0
k=0 βk q
Because L(q, θo ) and F (q, θo ) are co-prime by Assumption 1
and polynomials of order ml − 1 and mf , and α(q) and
β(q) are polynomials of orders at most mf − 1 and ml − 1,
(B) can only be satisfied if α(q) ≡ 0 ≡ β(q). Hence,
[−Qf∞ (ηo ) Ql∞ (ηo )] has full column rank.
Analogously for [−Qc∞ (ηo ) Qd∞ (ηo )], this matrix has full
column rank if and only if C(q, θo )α(q) = D(q, θo )β(q) is
satisfiedP
only for α(q) ≡ 0 ≡ β(q),
here we have
Pmdwhere
mc −1
−1
α(q) = k=0
αk q −k and β(q) = k=0
βk q −k . This is the
case, as C(q, θo ) and D(q, θo ) are co-prime and polynomials of
higher order than α(q) and β(q). Hence, [−Qf∞ (ηo ) Ql∞ (ηo )]
and [−Qc∞ (ηo ) Qd∞ (ηo )] are full column rank, implying that
Q∞ (ηo ) has a trivial right null-space and M (ηo ) is invertible.

Returning to (B), we may now write
LS
θ̂N
− θo ≤ C η̂N − η̂on(N ) , ∀N > N̄ .
→ 0, as N → ∞, w.p.1.

(59)

Moreover, using (B), we can re-write (B) as
LS
θ̂N
− θo ≤ C

η̂N − η̄ n(N ) + η̄ n(N ) − ηon(N )
n(N )


.

From Proposition 1, we have ||η̄ n(N ) − ηo
|| ≤ Cd(N ),
which thus approaches zero faster than ||η̂N − η̄ n(N ) ||, whose
LS
decay rate is according to (2). For the decay rate of ||θ̂N
−θo ||,
it suffices then to take the rate of the slowest-decaying term,
which is given by (1), as we wanted to show.
A PPENDIX C
C ONSISTENCY OF S TEP 3
The main purpose of this appendix is to prove Theorem 2.
However, before we do so, we introduce some results regarding the norm of some vectors and matrices.
n
• kRN k is bounded for all n and sufficiently large N , w.p.1
n
Let RN
be defined as in (II-D). Then, from Lemma 4.2 in [41],
we have that there exists N̄ such that, w.p.1,

Finally, we have the necessary results to prove Theorem 1.
Proof of Theorem 1: We start by using (III) to write

−1 ⊤
LS
θ̂N
−θo = Q⊤
Qn(η̂N)η̂N − θo
n(η̂N)Qn (η̂N)
n
kRN
k ≤ C, ∀n, ∀N > N̄ .
(60)
 ⊤
− 1 ⊤
= Qn(η̂N)Qn (η̂N) Qn(η̂N) [η̂N − Qn (η̂N)θo ]
−1
• ||Tn (θo )|| is bounded for all n
 ⊤
− 1 ⊤
n(N )
= Qn(η̂N)Qn (η̂N) Qn(η̂N)Tn (θo )[η̂N −ηo
]
We observe that, with Tn (θ) is given by (III), the inverse of
(56)
Tn (θo ) is given by
where the last equality follows from (III). If n were fixed,


consistency would follow if η̂N − ηon would approach zero as
Tnc (θ)−1
0
−1 (θ) =
,
(61)
T
n
N → ∞, provided the inverse of Q⊤
Tnf (θ)−1 Tnl (θ)Tnc (θ)−1 Tcf (θ)−1
n (η̂N )Qn (η̂N ) existed for
sufficiently large N . However, n = n(N ) increases according
to Assumption 4. This implies that the dimensions of the evaluated at the true parameters θo . Also, Tnf (θo )−1 , Tnc (θo )−1 ,
n(N )
, and of the matrices Qn (η̂N ) (number and Tnl (θo ) are sub-matrices of T[1/Fo (q)], T[1/Co (q)], and
vectors η̂N and ηo
of rows) and Tn (θo ) (number of rows and columns), become T[Lo (q)], respectively, where T[X(q)] is defined by (B). Then,
arbitrarily large. Therefore, extra requirements are necessary.
Tn−1 (θo ) ≤ Tnf (θo )−1 + Tnc (θo )−1
In particular, we use (B) to write
+ Tnf (θo )−1 Tnl (θo ) Tnc (θo )−1
n(N )
LS
(η̂
)T
(θ
)(η̂
−
η
)
θ̂N
− θo = M −1 (η̂N )Q⊤
N
n o
N
n
o
≤ kT[1/Fo (q)]k + kT[1/Co (q)]k
≤ M −1 (η̂N ) Qn (η̂N ) Tn (θo ) η̂N − ηon(N ) , (57)
+ kT[1/Fo (q)]k kT[Lo (q)]k kT[1/Co (q)]k
where M (η̂N ) := Q⊤
(η̂
)Q
(η̂
).
Consistency
is
achieved
≤C, ∀n,
(62)
N
n N
n
if the last factor on the right side of the inequality in (B)
approaches zero, as N → ∞, w.p.1, and the remaining factors where the last inequality follows from (B) and the fact that
are bounded for sufficiently large N , w.p.1. This can be shown 1/Fo (q), 1/Co (q), and Lo (q) are stable transfer functions.
LS
)|| is bounded for all n and sufficiently large N
using (B), (B), and (B), but we need additionally that M (η̂N ) • ||Tn−1 (θ̂N
LS
LS
is invertible for sufficiently large N , w.p.1.
)||, where Tn−1 (θ̂N
) is given by (C)
Consider the term ||Tn−1 (θ̂N
LS
With this purpose, we write
evaluated at θ̂N . We have that, proceeding as in (C),
||M (η̂N ) − M (ηon(N ) )||
⊤ n(N )
= ||Q⊤
)Qn (ηon(N ) )||
n (η̂N )Qn (η̂N ) − Qn (ηo

LS
LS
LS
) ≤ T[1/C(q, θ̂N
)] + T[1/F (q, θ̂N
)]
Tn−1 (θ̂N

Using (B), (B), (B), and Proposition 5, and because
M (ηon ) → M (ηo ) as n → ∞, we have that
M (η̂N ) → M (ηo ), as N → ∞, w.p.1.

(58)

As M (ηo ) is invertible (Lemma 1), by (B) and because the map
from the entries of a matrix to its eigenvalues is continuous,
there is N̄ such that M (η̂N ) is invertible for all N > N̄ , w.p.1.

LS
+ T[1/C(q, θ̂N
)]

LS
T[L(q, θ̂N
)]

LS
T[1/F (q, θ̂N
)]

LS
LS
for all n. This will be bounded if F (q, θ̂N
) and C(q, θ̂N
) have
all poles strictly inside the unit circle. From Theorem 1 and
stability of the true system by Assumption 1, we conclude that
LS
LS
there exists N̄ such that F (q, θ̂N
) and C(q, θ̂N
) have all roots
strictly inside the unit circle for all N > N̄ . Thus, we have
LS
) ≤ C, ∀n, ∀N > N̄ , w.p.1.
Tn−1 (θ̂N

(63)

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

13

LS
||Tn−1 (θ̂N
) − Tn−1 (θo )|| tends to zero, as N tends to infinity,
w.p.1
LS
For the term ||Tn−1 (θ̂N
) − Tn−1 (θo )||, with n = n(N ), we have

•

LS
) − Tn−1 (θo ) ≤
Tn−1 (θ̂N
(64)
LS
LS
≤ Tn−1 (θ̂N
) Tn (θ̂N
) − Tn (θo ) Tn−1 (θo ) .
p
Because ||X|| ≤
||X||1 ||X||∞ , where X is an arbitrary
matrix, we have that
Pmf +ml +mc k,LS
LS
Tn (θ̂N
) − Tn (θo ) ≤ k=1
|θ̂N − θok |
LS
− θo ,
≤ C θ̂N

with superscript k denoting the k th element of the vector, we
can use Theorem 1 to show that
r


log N
LS
n(N )
Tn (θ̂N )− Tn(θo ) = O
1 + d(N ) . (65)
N
From (IV) and (IV),
r

log N
1 + d(N ) → 0, as N → ∞,
n2 (N )
N

Together with (C), (C), and (C), this implies that
(66)

The following two lemmas are useful for the invertibility of
the weighted least squares problem (III).
Lemma 2. Let Assumption 1 hold and
n→∞

(67)

where W̄n (θo ) is given by (III), and Qn (ηon ) is defined by (III)
at the true parameters ηon . Then, M̄ (ηo , θo ) is invertible.
Proof. Using (II-D) and (III), we re-write (2) as
M̄ (ηo , θo ) =
lim Q⊤ (ηon )Tn−⊤ (θo )Ē
n→∞ n



ϕnt (ϕnt )⊤

 −1
Tn (θo )Qn (ηon ). (68)

Re-writing ϕnt , defined in (II-D), as
 

 
−Γn Go (q) −Γn Ho (q) ut
−Γn yt
=
,
ϕnt =
Γn
0
Γn u t
et
we can then write
Z π
 n n ⊤
1
Ē ϕt (ϕt ) =
Λn (eiω )Φz Λ∗n (eiω )dω,
2π −π
where

0

where the argument q of the polynomials was dropped for
notational simplicity. In turn, we can also write
Q⊤
(η n )T −⊤ (θo )Λn =
n o n
 
Bo
⊤
−Tn,m
Γn

f Fo

Ao
⊤
Γn
T




  n,ml Fo
 ⊤ Ao
L o Ao
⊤
Γn
 Tn,mc Co Γn Go−Tn,m
c Fo Co



 
Lo
1
⊤
⊤
−Tn,md Co Γn Go+Tn,md Fo Co Γn

n→∞

LS
Tn (θ̂N
) − Tn (θo ) → 0, as N → ∞, w.p.1.

n
n
M̄ (ηo , θo ) := lim Q⊤
n (ηo )W̄n (θo )Qn (ηo ),

 
Bo
⊤
−Tn,m
f

 Fo 


Ao
⊤
0
Tn,m


l
Fo
⊤ n
−
⊤
,





Qn (ηo )Tn (θo ) = 
L o Ao 
Ao
⊤
⊤
−Tn,mc Fo Co 
−Tn,mc Co


 

Lo
1
⊤
⊤
Tn,m
Tn,m
Co
Fo Co
d
d




0




.
 

Ao
⊤
Tn,m
Γ
H

n
o
c Co

 
1
⊤
−Tn,md Co Γn Ho
0

It is possible
to observe that, for some polynomial
P∞
⊤
X(q) = k=0 xk q −k , lim Tn,m
(X(q))Γn = X(q)Γm . Then,
n→∞
using also (II-D), we have lim Q⊤ (η n )T −⊤ (θ )Λ = Ω,

and thus

LS
) − Tn−1 (θo ) → 0, as N → ∞, w.p.1.
Tn−1 (θ̂N

Moreover,




−Γn Go (q) −Γn Ho (q)
Λn (q) =
.
Γn
0

Then, we can re-write (C) as
Z π
1
lim Q⊤ (η n )T −⊤ (θo )Λn (eiω )
M̄ (ηo , θo ) =
2π −π n→∞ n o n
· Φz Λ∗n (eiω )Tn−1 (θo )Qn (ηon )dω. (69)

n

o

n

o

n

where Ω is given by (II-C). This allows us to re-write (C)
as
Z π
1
ΩΦz Ω∗ dω = MCR ,
(70)
M̄ (ηo , θo ) =
2π −π
which is invertible because the CR bound exists for an
informative experiment [1].

Lemma 3. Let Assumptions 1, 2, 3, and 4 hold.
LS
LS
Also, let M (η̂N , θ̂N
) := Q⊤
n (η̂N )Wn (θ̂N )Qn (η̂N ), where
n(N )
η̂N := η̂N
is defined by (IV), Qn (η̂N ) is defined by (III)
LS
evaluated at the estimated parameters η̂N , and Wn (θ̂N
) is
defined by (III). Then,
LS
M (η̂N , θ̂N
) → M̄ (ηo , θo ), as N → ∞, w.p.1.

(71)

Proof. For the purpose of showing (3), we will show that
LS
n(N )
M (η̂N , θ̂N
) − Q⊤
)W̄n (θo )Qn (ηon(N ) )
n (ηo
→ 0, as N → ∞, w.p.1.

(72)

We apply Proposition 5, whose conditions can be verified
using (B), (B), (C), and (C), but additionally we need to verify
LS
that ||W̄n (θo )|| is bounded and ||Wn (θ̂N
) − W̄n (θo )|| tends to
zero. For the first, we have that
||W̄n (θo )|| ≤ ||Tn−1 (θo )||2 ||R̄n || ≤ C,

(73)

LS
following from (C) and (C). For the second, with Wn (θ̂N
)
given by (III) and W̄n (θo ) by (III), conditions in Proposition 5
are satisfied using (C), (C), (C) and Proposition 3, from where
it follows that
LS
||Wn (θ̂N
) − W̄n (θo )|| → 0, as N → ∞, w.p.1.

(74)

Having shown (C) and (C), the assumptions of Proposition 5
are verified, from which (C) follows and implies (3).
We now have the necessary results to prove Theorem 2.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

14

Proof of Theorem 2: Similarly to (B), we write
WLS
θ̂N
− θo
LS
LS
n(N )
= M −1 (η̂N , θ̂N
)Q⊤
),
n (η̂N )Wn (θ̂N )Tn (θo )(η̂N − ηo
(75)

and analyze
WLS
LS
θ̂N
− θo ≤ M −1 (η̂N , θ̂N
)

·

LS
Wn (θ̂N
)

Qn (η̂N )
Tn (θo )

η̂N − ηon(N ) .

and δ̂N is a stochastic vector of compatible dimensions. The
dimensions may increase to infinity as function of N , except for
the number of rows of ÂN , which is fixed. We assume that there
is N̄ such that ||ÂN || < C for all N > N̄ , there is B̄ such that
||B̂N − B̄|| → 0√as N → ∞ w.p.1, and ||δ̂N || → 0 as N → ∞
w.p.1. Then,
√ if N ||B̂N − B̄|| ||δ̂N || → 0 as N → ∞ w.p.1,
x̂N and N ÂN B̄ δ̂N have the same asymptotic distribution
and covariance.
Proof. We begin by writing
√
√
x̂N = N ÂN B̄ δ̂N + N ÂN (B̂N − B̄)δ̂N .

LS
(77)
From Lemma 3, M (η̂N , θ̂N
) converges to M̄ (ηo , θo ), which
is invertible from Lemma 2. Hence, because the map from Although some of the matrix and vector dimensions may
the entries of the matrix to its eigenvalues is continuous, increase to infinity with N , the number of rows of ÂN is fixed,
LS
M (η̂N , θ̂N
) is invertible for sufficiently large N , and therefore which makes x̂N finite dimensional, to which [51, Lemma B.4]
its norm is bounded, as it is a matrix of fixed dimensions. may be applied. Then, it is a consequence of this lemma that
√
Also, from (B), kQn (η̂N )k is bounded for sufficiently large x̂N and N ÂN B̄ δ̂N have the same asymptotic distribution
N . Moreover, we have that, making explicit that n = n(N ),
and covariance if the second term on the right side of (D)
n(N )
−
1
LS 2
LS
tends to zero with probability one. By assumption, we have
RN
.
Wn(N ) (θ̂N ) ≤ Tn(N ) (θ̂N )
√
√
|| N ÂN (B̂N − B̄)δ̂N || ≤ N ||ÂN || ||B̂N − B̄|| ||δ̂N ||
Then, from (C) and (C), we have
→ 0, as N → ∞, w.p.1,
LS
Wn(N ) (θ̂N
) ≤ C, ∀N > N̄ .

which completes the proof.

Finally, using also (B), (B), and (B), we conclude that
WLS
θ̂N
− θo → 0, as N → ∞, w.p.1.

A PPENDIX D
A SYMPTOTIC D ISTRIBUTION AND C OVARIANCE

We now have the necessary results to prove Theorem 3.
Proof of Theorem 3: We start by re-writing (D) as
√
WLS
LS
LS
N (θ̂N
− θo ) = M −1 (η̂N , θ̂N
)x(η̂N , θ̂N
),
where

OF

S TEP 3

The purpose of this appendix is to prove Theorem 3:
asymptotic distribution and covariance of
√
√
WLS
LS
N (θ̂N
− θo ) = N Υn (η̂N , θ̂N
)(η̂N − ηon(N ) ), (76)
which is re-written from (C), where
LS
LS
−1
Υn (η̂N , θ̂N
) = [Q⊤
n (η̂N )Wn (θ̂N )Qn (η̂N )]
⊤
LS
· Qn (η̂N )Wn (θ̂N )Tn (θo ).
LS
If Υn (η̂N , θ̂N
) were of fixed dimensions, the standard idea
LS
would be to show that Υn (η̂N , θ̂N
) converges w.p.1 to a
LS
being
deterministic matrix, as consequence of η̂N and θ̂N
consistent estimates of ηo and θo , respectively. Then, for
computing the asymptotic distribution and covariance of (D),
one√can consider the asymptotic distribution and covariance
n(N )
LS
of N (η̂N − ηo
) while Υn (η̂N , θ̂N
) can be replaced by
the deterministic matrix it converges to. This standard result
follows from [51, Lemma B.4], but it is not applicable here
n(N )
LS
because the dimensions of Υn (η̂N , θ̂N
) and η̂N − ηo
are
not fixed. In this scenario, Proposition 4 must be used instead.
However, (D) is not ready to be used with Proposition 4,
because it requires η̂N − ηon to be pre-multiplied by a deterministic matrix. The key idea of proving Theorem 3 is to show
that (D) has the same asymptotic distribution and covariance
n(N )
as an expression of the form Ῡn [η̂N − ηo
], where Ῡn
is a deterministic matrix, and then apply Proposition 4. The
following result will be useful for this purpose.
√
N ÂN B̂N δ̂N be a finiteProposition 6. Let x̂N =
dimensional vector, where ÂN and B̂N are stochastic matrices

LS
LS
M (η̂N , θ̂N
) = Q⊤
n (η̂N )Wn (θ̂N )Qn (η̂N ),
√
LS
n(N )
LS
).
x(η̂N , θ̂N
) = N Q⊤
n (η̂N )Wn (θ̂N )Tn (θo )(η̂N −ηo

LS
LS
Both M (η̂N , θ̂N
) and x(η̂N , θ̂N
) are of fixed dimension, and
we have from (C) and (3) that
−1 , as N → ∞, w.p.1.
LS
M −1 (η̂N , θ̂N
) → MCR

Then, if we assume that
LS
x(η̂N , θ̂N
) ∼ AsN (0, P ),

we have that, from [51, Lemma B.4],


√
−
1
−
1
WLS
N (θ̂N − θo ) ∼ AsN 0, MCR P MCR .

(78)

(79)

We will proceed to show that (D) is verified with

n
n
2
P = σo2 lim Q⊤
n (ηo )W̄n (θo )Qn (ηo ) = σo MCR ,
n→∞

(80)

where the second equality follows directly from (2) and (C).
We now proceed to show the first equality.
In the following arguments, we will apply Proposition 6
LS
repeatedly to x(η̂N , θ̂N
). This required the boundedness of
LS
some matrices; however, because all the matrices in x(η̂N , θ̂N
)
have been shown to be bounded for sufficiently large N w.p.1,
for readability we will refrain from referring to this every time
Proposition 6 is applied.
Because it is more convenient to work with η̄ n(N ) than
n(N )
LS
ηo
, we start by re-writing x(η̂N , θ̂N
) as
√
⊤
LS
LS
x(η̂N , θ̂N ) = N Qn (η̂N )Wn (θ̂N )Tn√(θo )(η̂N − η̄ n(N ) )
LS
n(N )
+ Q⊤
− ηon(N ) ).
n (η̂N )Wn (θ̂N )Tn (θo ) N (η̄

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

Using Proposition 1, we have, for sufficiently large N w.p.1,
√
LS
n(N )
Q⊤
− ηon(N ) ) ≤
n (η̂N )Wn (θ̂N )Tn (θo ) N (η̄
√
≤ C N d(N ) → 0, as N → ∞.
Using an identical argument to Proposition 6, we have that
LS
x(η̂N , θ̂N
) and
√
LS
n(N )
N Q⊤
)
(81)
n (η̂N )Wn (θ̂N )Tn (θo )(η̂N − η̄
have the same asymptotic distribution and covariance, so we
will analyze (D) instead.
LS
Expanding Wn (θ̂N
) in (D), we obtain
√
⊤
LS
NQ
)Tn (θo )(η̂N − η̄ n(N ) )
√n (η̂N⊤)Wn (θ̂N
−⊤ LS
n −1 LS
= N Qn(η̂N )Tn (θ̂N )RN
Tn (θ̂N )Tn (θo )(η̂N − η̄ n(N ) ).
(82)
Using Proposition 6, we conclude that (D) and
√
−⊤ LS
n −1
N Q⊤
(θ̂N )RN
Tn (θo )Tn (θo )(η̂N − η̄ n(N ) )
n(η̂N )Tn √
⊤
LS
n
)RN
(η̂N − η̄ n(N ) )
= N Qn (η̂N )Tn−⊤ (θ̂N
(83)
have the same asymptotic properties if
√
LS
N ||Tn−1 (θ̂N
) − Tn−1 (θo )|| ||η̂N − η̄ n(N ) || → 0,
as N → ∞, w.p.1.
Using (C), (C) and (C) to write
√
LS
−1
n(N )
N ||Tn−1 (θ̂N
||
o )|| ||η̂N − η̄
√ ) − Tn (θ
LS
≤ C N ||Tn (θ̂N ) − Tn (θo )|| ||η̂N − η̄ n(N ) ||,
we have from (C) and Proposition 2 that
√
LS
N Tn (θo ) − Tn (θ̂N
) η̂N − η̄ n(N ) =


2
n(N ) log N
√
=O
, (84)
1 + d(N )
N

where

n(N ) log N
√
=
N



n3+δ (N )
N

1
 3+δ

log N

→ 0, as N → ∞,
1+δ
N 2(3+δ)
due to Condition D2 in Assumption 4. This implies that (D)
LS
and (D), and in turn x(η̂N , θ̂N
), have the same asymptotic
distribution and covariance. Repeating this procedure, it can
be shown that these, in turn, have the same asymptotic
distribution and covariance as
√
−⊤
n
n(N )
N Q⊤
).
(85)
n (η̂N )Tn (θo )RN (η̂N − η̄
There are two stochastic matrices left in (D), which we
need to replace by deterministic matrices that do not affect
the asymptotic properties. Using Proposition 3,
√
n
N RN
− R̄n η̂N − η̄ n(N ) =
 3/2

n (N ) log N
√
1 + d(N )
=O 2
N
r
r


n2 (N ) log N n3 (N )
1 + d(N ) ,
+C
N
N

15

where the first term tends to zero by applying Condition D2
in Assumption 4 to
n3/2 (N ) log N
√
=
N



n4+δ (N )
N

3
 2(4+δ)

log N
1+δ

N 2(4+δ)
→ 0, as N → ∞,

and the second because of Condition D2 in Assumption 4,
and (IV). Then, from Proposition 6, we have that (D) and
√
−⊤
n
n(N )
N Q⊤
),
(86)
n (η̂N )Tn (θo )R̄ (η̂N − η̄
LS
and in turn x(η̂N , θ̂N
), have the same asymptotic distribution
and covariance, so we will analyze (D).
Applying again Proposition 6, we have that (D) and
√
N Q⊤ (η̄ n(N ) )T −⊤ (θ )R̄n (η̂ − η̄ n(N ) )
(87)
n

n

o

N

have the same asymptotic properties, since
√
N Qn (η̂N ) − Qn (η̄ n(N ) ) η̂N − η̄ n(N )
√
2
≤ C N η̂N − η̄ n(N )


2
n(N ) log N
√
=O
1 + d(N )
N
by using (B) and Proposition 2, which tends to zero as
N → ∞, identically to (D).
In (D), the matrix multiplying η̂N − η̄ n(N ) is finally detern(N )
ministic, but it will be more convenient to work with Q(ηo
).
With this purpose, Proposition 1 can be used to show that (D)
and
√
n(N )
N Q⊤
)Tn−⊤ (θo )R̄n (η̂N − η̄ n(N ) )
(88)
n (ηo
have the same asymptotic properties, as
√
N Qn (η̄ n(N ) ) − Qn (ηon ) η̂N − η̄ n(N )
√
≤ C N η̄ n(N ) − ηon(N ) η̂N − η̄ n(N )

r
√
n(N ) log N
1 + d(N ) N d(N ) ,
=O
N

LS
which tends to zero due to (IV) and (IV). Thus, x(η̂N , θ̂N
)
and (D) have the same asymptotic distribution and covariance,
so we will analyze (D) instead.
n
−⊤
n
Let Υn := Q⊤
n (ηo )Tn (θo )R̄ . Then, using Proposition 4,
√
N Υn (η̂N − η̄ n(N ) ) ∼ AsN (0, P ),

where P is given by (D). Finally, using (C), (D), and (D):
√
−1 ).
WLS
N (θ̂N
− θo ) ∼ AsN (0, σo2 MCR
R EFERENCES
[1] L. Ljung System Identification. Theory for the User. Prentice-Hall, 1999.
[2] K. Åström and T. Bohlin. Numerical identification of linear dynamic
systems from normal operating records. In IFAC Symposium on SelfAdaptive Systems, Teddington, United Kingdom, 1965.
[3] T. Söderström and P. Stoica. Instrumental variable methods for system
identification. Springer-Verlag, 1983.
[4] P. Stoica and T. Söderström. Optimal instrumental variable estimation
and approximate implementations. IEEE Transactions on Automatic
Control, vol. 28:757–772, July 1983.
[5] P. Young. The refined instrumental variable method: Unified estimation
of discrete and continuous-time transfer function models. Journal
Européen des Systèmes Automatisés, 42:149–179, 2008.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

[6] M. Gilson and P. van den Hof. Instrumental variable methods for closedloop system identification. Automatica, 41(2):241–249, 2005.
[7] S. Y. Kung. A new identification and model reduction algorithm via
singular value decomposition. In Asilomar Conference on Circuits,
Systems and Computers, pages 705–714, Pacific Grove, USA, 1978.
[8] P. van Overschee and B. de Moor. Subspace identification for linear systems: theory—implementation—applications. Kluwer Academic, Boston,
1996.
[9] M. Verhaegen. Application of a subspace model identification technique to identify LTI systems operating in closed-loop. Automatica,
29(4):1027–1040, 1993.
[10] S. Qin and L. Ljung. Closed-loop subspace identification with innovation
estimation. In 13th IFAC Symposium on System Identification, pages
861–866, Rotterdam, The Netherlands, 2003.
[11] M. Jansson. Subspace identification and ARX modeling. In 13th IFAC
Symposium on System Identification, Rotterdam, The Netherlands, 2003.
[12] A. Chiuso and G. Picci. Consistency analysis of some closed-loop
subspace identification methods. Automatica, 41(3):377–391, 2005.
[13] M. Jansson and B. Wahlberg. A linear regression approach to state-space
subspace system identification. Signal Processing, 52:103–129, 1996.
[14] D. Bauer. Asymptotic properties of subspace estimators. Automatica,
41(3):359 – 376, 2005.
[15] C. Sanathanan and J. Koerner. Transfer function synthesis as a ratio
of two complex polynomials. IEEE Transactions on Automatic Control,
8(1):56–58, 1963.
[16] A. G. Evans and R. Fischl. Optimal least squares time-domain synthesis
of recursive digital filters. IEEE Transactions on Audio and Electroacoustics, 21(1):61–65, 1973.
[17] Y. Bresler and A. Macovski. Exact maximum likelihood parameter estimation of superimposed exponential signals in noise. IEEE Transactions
on Signal Processing, 34(5):1081–1089, 1986.
[18] A. K. Shaw. Optimal identification of discrete-time systems from
impulse response data. IEEE Trans. on Signal Processing, 42(1):113–
120, 1994.
[19] A. K. Shaw, P. Misra, and R. Kumaresan. Identification of a class of
multivariable systems from impulse response data: Theory and computational algorithm. Circuits, Systems and Signal Processing, 13(6):759–
782, 1994.
[20] P. Lemmerling, L. Vanhamme, S. van Huffel, and B. de Moor. IQMLlike algorithms for solving structured total least squares problems: a
unified view. Signal Processing, 81:1935–1945, 2001.
[21] K. Steiglitz and L. E. McBride. A technique for the identification of
linear systems. IEEE Trans. on Automatic Control, 10:461–464, 1965.
[22] J. H. McClellan and D. Lee. Exact equivalence of the SteiglitzMcBride iteration and IQML. IEEE Transactions on Signal Processing,
39(2):509–512, 1991.
[23] P. Stoica and T. Söderström. The Steiglitz-McBride identification
algorithm revisited–convergence analysis and accuracy aspects. IEEE
Transactions on Automatic Control, 26(3):712–717, 1981.
[24] T. Söderström, P. Stoica, and B. Friedlander. An indirect prediction error
method for system identification. Automatica, 27(1):183–188, 1991.
[25] C. L. Byrnes, T. Georgiou, and A. Lindquist. A new approach to
spectral estimation: A tunable high-resolution spectral estimator. IEEE
Transactions on Signal Processing, 48(11):3189–3205, 2000.
[26] M. Zorzi. An interpretation of the dual problem of the three-like
approaches. Automatica, 62:87–92, 2015.
[27] E. J. Hannan. The asymptotic theory of linear time-series models.
Journal of Applied Probability, 10(1):130–145, 1973.
[28] J. Durbin. The fitting of time-series models. Revue de l’Institut
International de Statistique, pages 233–244, 1960.
[29] James Durbin. Efficient estimation of parameters in moving-average
models. Biometrika, 46(3/4):306–316, 1959.
[30] D. Mayne and F. Firoozan. Linear identification of ARMA processes.
Automatica, 18(4):461–466, 1982.
[31] E. J. Hannan and L. Kavalieris. Linear estimation of ARMA processes.
Automatica, 19(4):447–448, 1983.
[32] E. Hannan and L. Kavalieris. Multivariate linear time series models.
Advances in Applied Probability, 16(3):492–561, 1984.
[33] G. Reinsel, S. Basu, and S. Yap. Maximum likelihood estimators in the
multivariate autoregressive moving-average model from a generalized
least squares viewpoint. J. Time Series Anal., 13(2):133–145, 1992.
[34] J. M. Dufour and T. Jouini. Asymptotic distributions for quasi-efficient
estimators in echelon VARMA models. Computational Statistics & Data
Analysis, 73:69–86, 2014.
[35] B. Wahlberg. Model reduction of high-order estimated models: the
asymptotic ML approach. Int. Journal of Control, 49(1):169–192, 1989.

16

[36] Y. Zhu. Multivariable System Identification for Process Control. Pergamon, 2001.
[37] Y. Zhu and H. Hjalmarsson. The Box-Jenkins Steiglitz-McBride algorithm. Automatica, 65:170–182, 2016.
[38] N. Everitt, M. Galrinho, and H. Hjalmarsson. Open-loop asymptotically
efficient model reduction with the Steiglitz–McBride method. Automatica, 89:221–234, 2018.
[39] M. Galrinho, C. R. Rojas, and H. Hjalmarsson. A weighted least-squares
method for parameter estimation in structured models. In IEEE Conf.
on Decision and Control, pages 3322–3327, Los Angeles, USA, 2014.
[40] U. Forssell and L. Ljung. Closed-loop identification revisited. Automatica, 35:1215–1241, 1999.
[41] L. Ljung and B. Wahlberg. Asymptotic properties of the least-squares
method for estimating transfer functions and disturbance spectra. Advances in Applied Probabilities, 24:412–440, 1992.
[42] C. Gourieroux and A. Monfort. Statistics and econometric models,
volume 2. Cambridge University Press, 1995.
[43] T. Kailath, A. H. Sayed, and B. Hassibi. Linear Estimation. PrenticeHall, 2000.
[44] M. Galrinho, C. R. Rojas, and H. Hjalmarsson. Estimating models with
high-order noise dynamics using semi-parametric Weighted Null-Space
Fitting. to be submitted (arXiv 1708.03947).
[45] M. Galrinho, C. R. Rojas, and H. Hjalmarsson. On estimating initial
conditions in unstructured models. In IEEE Conference on Decision and
Control, pages 2725–2730, Osaka, Japan, 2015.
[46] W. Larimore. System identification, reduced-order filtering and modeling
via canonical variate analysis. In American Control Conference, 1983.
[47] K. Peternell, W. Scherrer, and M. Deistler. Statistical analysis of novel
subspace identification methods. Signal Processing, 52:161–177, 1996.
[48] P. van Overschee and B. de Moor. N4SID: Subspace algorithms for the
identification of combined deterministic-stochastic systems. Automatica,
30:75–93, 1994.
[49] P. Stoica and M. Jansson. MIMO system identification: state-space
and subspace approximations versus transfer function and instrumental
variables. IEEE Trans. on Signal Processing, 48(11):3087–3099, 2002.
[50] V. Peller. Hankel operators and their applications. Springer, 2003.
[51] T. Söderström and P. Stoica. System Identification. Prentice Hall, 1989.

PLACE
PHOTO
HERE

Miguel Galrinho was born in 1988. He received his
M.S. degree in aerospace engineering in 2013 from
Delft University of Technology, The Netherlands,
and the Licentiate degree in electrical engineering
in 2016 from KTH Royal Institute of Technology,
Stockholm, Sweden.
He is currently a PhD student at KTH, with the
Department of Automatic Control, School of Electrical Engineering, under supervision of Professor
Håkan Hjalmarsson. His research is on least-squares
methods for identification of structured models.

PLACE
PHOTO
HERE

Cristian R. Rojas (M’13) was born in 1980. He
received the M.S. degree in electronics engineering
from the Universidad Técnica Federico Santa Marı́a,
Valparaı́so, Chile, in 2004, and the Ph.D. degree in
electrical engineering at The University of Newcastle, NSW, Australia, in 2008.
Since October 2008, he has been with the Royal
Institute of Technology, Stockholm, Sweden, where
he is currently Associate Professor at the Department
of Automatic Control, School of Electrical Engineering. His research interests lie in system identification

and signal processing.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

Håkan Hjalmarsson (M’98–SM’11–F’13) was
born in 1962. He received the M.S. degree in electrical engineering in 1988, and the Licentiate and
Ph.D. degrees in automatic control in 1990 and
PLACE
1993, respectively, all from Linköping University,
PHOTO
Linköping, Sweden.
HERE
He has held visiting research positions at California Institute of Technology, Louvain University,
and at the University of Newcastle, Australia. His research interests include system identification, signal
processing, control and estimation in communication
networks and automated tuning of controllers.
Dr. Hjalmarsson has served as an Associate Editor for Automatica (1996–
2001) and for the IEEE Transactions on Automatic Control (2005–2007),
and has been Guest Editor for the European Journal of Control and Control
Engineering Practice. He is a Professor at the School of Electrical Engineering,
KTH, Stockholm, Sweden. He is a Chair of the IFAC Coordinating Committee
CC1 Systems and Signals. In 2001, he received the KTH award for outstanding
contribution to undergraduate education. He is co-recipient of the European
Research Council advanced grant.

17

