DietmarBauer
Departmentof ElectricalEngineering
Linkoping University, SE-58183 Linkoping,Sweden
WWW: http://www.control.isy.l iu.s e
Email: Dietmar.Bauer@tuwien.ac.at
June 6,2000
REG
LERTEKNIK
AUTO
MATIC CONTR
OL
LINKÖPING
Report no.: LiTH-ISY-R-2265
Submitted toCDC'2000, Sydney, Dec. 2000
TechnicalreportsfromtheAutomaticControlgroupinLinkopingareavailablebyanonymousftp
DietmarBauer
Institutef. Econometrics, OperationsResearch and System Theory,
TUWien, Argentinierstr. 8,A-1040 Wien
June 6,2000
Abstract
Inmoderndata analysisoften the rststepis to performsomedata preprocessing, e.g.
detrendingor eliminationofperiodiccomponentsofknownperiodlength. Thisis normally
done using leastsquares regression. Onlyafterwards black box models are estimatedusing
eitherpseudo-maximum-likelihoodmethods,predictionerrormethodsorsubspacealgorithms.
Inthispaperitisshown,thatfor subspacemethodsthisisessentiallythesameasincluding
the corresponding input variables, e.g. a constant or a trend or aperiodic component, as
additionalinputvariables. Hereessentiallymeans,thattheestimatesonlydierthroughthe
choiceofinitialvalues.
Keywords: subspacealgorithms,linearsystems,estimation,identication
1 Introduction
Ithasbecomestandardinmoderntimeseriesanalysistoperformsomeform ofpreprocessingon
thedatapriortoidentication. Howeveroftentheeectsofthispreprocessingarenotdealt with
in the identication phase. As an exampleconsider the analysis of so called subspacemethods:
Manyalgorithmshavebeenproposedbydierentauthors,whichareallsubsumedunderthename
'subspacealgorithms',asthey showcertain similarities. Toname just themost popular wecite
CCA(Larimore, 1983), MOESP(Verhaegen, 1994)and N4SID(VanOverscheeand DeMoor,1994).
Thepropertiesof thesealgorithms havebeen analysedin anumberof papers(forreferencessee
e.g.Baueretal.,1999). Howeverallcitedreferencesrefertothecase,thatnopreprocessingprior
toidenticationtakesplace. Alsonormallysomepersistencyconditionsontheinputareimposed,
which excludese.g. the constant as an inputvariable. Therefore essentiallyit is assumed, e.g.
that theonly sourcefor anonzero mean of the output is due to theltered input. Whenthis
doesnotholdandnopreprocessingisperformed,abiaswillresult. Ifpreprocessingisperformed
the properties of thedata change, which has to be takeninto account forthe estimation of the
uncertainties, i.e. for the estimation of the asymptotic variances of the model obtained by the
identication step. This in many cases will be straightforward, however it seems to be more
convenient to perform both thepreprocessing and theactual identicationin one step,in order
tounify thetreatmentandtosimplifythecalculationoftheasymptoticvariance.
Theaimofthispaperistoshow,thatthepreprocessingandtheestimationphaseinfactcanbe
uniedbyincludingtheconstantorsimilartermsasinputvariablesandusegeneralisedinverses
i.e. regularisationtechniquesintherespectiveregressions. Infact,itwillbeshown,thatthisleads
toessentiallythesameresultsinthesense,thatthereisnodierenceintheattainedaccuracy,i.e.
that theasymptoticdistributionof theestimated systemareidentical. Thisshows,that someof
thepreprocessingmightbeincludedintothesubspacealgorithmwithoutanyproblem,and thus
canalsobeanalysedalongthesamelines.
Thiswork has been done, while the author was holdinga post doc positionat the Departmentfor
Auto-maticControl,UniversityofLinkoping,Sweden. Thenancialsupportbythe EUTMRproject'SI'isgratefully
arestated. Section3thengivesashortdescriptionofthemainstepsof theprocedure. Section4
then states the main results of the paper. In section 5 a numerical example is presented and
section6concludes.
2 Model set and assumptions
Inthispaperwedealwithdiscretetime,nitedimensional,timeinvariant,statespacesystemsof
theform x t+1 = Ax t +Bz t +K" t y t = Cx t +Dz t +" t (1) Herey t
denotesthesdimensionalobservedoutput,u
t
denotesthemdimensionalobservedinput,
"
t
the s dimensional white noise. A 2 R nn ;B 2 R nm ;C 2 R s n ;D 2 R s m ;K 2 R ns are
parameter matrices. The system usually will be described as (A;B;C ;D;K). We will assume
throughout, that "
t
is i.i.d. with zero mean and variance matrix > 0, having nite fourth
moments. Throughout wewillassumethat thesystemisstable, i.e. j
max
(A)j<1holds,where
max
denotes an eigenvalueof maximummodulus. It will also be assumed, that the systemis
strictlyminimum-phase,i.e. thatj
max
(A KC)j<1.
Correspondingtotheinputwewillassume,that theinput canbepartitionedinto twoparts:
u t 2 R mi and p t 2 R m mi . Here u t
accounts for the identication input, whereas p
t
denotes
theadditionalinputsdue to thepreprocessing. These additionalinputswill berestrictedto the
followingchoices:
aconstantterm,i.e. p
t =1;8t
aperiodiccomponent,i.e. p
t;1
=sin(!t+)forknown!2( ;]and. Inthiscasealso
thelaggedvariablep
t;2 =p
t 1;1
hastobeincluded
atimetrend,i.e. p
t =t;8t
These choices include the typicalpreprocessing likedetrendingand eliminating periodic
compo-nentsof known periodicity. Naturally,all these termscould be used in combination. Especially
the inclusion of the trend makes the inclusion of a constant necessary in order to account for
unknowninitial eects. Thekeyfeature ofthese inputsis that theyare persistent oforder one.
Thisisthereasonforincludingthelaggedtermfortheperiodiccomponents. Itwillbeclearfrom
thetext, why thisisneeded. Notethat in asimilar fashionalso morecomplicated preprocessing
canbe dealt with, aslong asit is done using regressiononto deterministic processes which are
persistentoforderone.
Note, thatthe model asstatedis notidentiable forthese inputs, which canbeseenfrom a
discussionoftheconstant: Assumethattheinputhasnonzeromean
u
andtheoutputhasmean
y . Thenweobtain y =D u + 1 X j=0 CA j B u
Thusitis clear,that theconstributionoftheconstantisequalto D
p + P 1 j=0 CA j B p , whereB p andD p
denotethecolumnsofB andD respectivelycorrespondingtotheconstant. Thisleadsto
alinearequationinB
m andD
m
showingthenonidentiability,whichwillcomplicatetheanalysis
oftheasymptoticdistribution.
Inthefollowingwewill always distinguishbetweenthese twokinds ofinputs mainly for the
reasonofstatingassumptionsfortheidenticationinputs. Thereasonisthedierentnatureofthe
inputs: All preprocessinginputsusedabovehavethecharacteristicofbeingperfectlypredictable
fromoneobservation,i.e. theyarepersistentlyexcitingofdegreeone. Fortheidenticationinputs
Subspaceprocedures havebeendescribedin anumberof papers. Thereforewerestrictourselves
to a short descriptiononly. Fordetails see the survey in (Bauer, 1998, Chapter 3). Let Y + t;f = [y 0 t ;;y 0 t+f 1 ] 0
denotethestackedvectorofoutputs,wheref isauserdenedinteger. Similarily
dene U + t;f and E + t;f
using the inputs and the noise respectively in the place of the outputs.
Furthermore letZ t;p =[y 0 t 1 ;u 0 t 1 ;;y 0 t p ;u 0 t p ] 0
, where pis auserdened integer. Then itis
easytoshowthefollowingequation:
Y + t;f =O f K p Z t;p +O f A p x t p +U f U + t;f +E f E + t;f : (2) Here A=(A KC), B=(B KD),O f =[C 0 ;A 0 C 0 ;:::;(A f 1 ) 0 ] 0 and K p =[[K ; B]; A [K ; B];:::; A p 1 [K ; B] FurtherE f
denotestheblockToeplitzmatrix,whosei-thblockrowisequalto
[CA i 2 K ;:::;CK ;I s ;0 s(f i)s ] where I s
is the s-dimensional identity and 0 ab
is a ab dimensional nullmatrix. Finally U
f
denotestheblockToeplitzmatrix,whosei-thcolumn isequalto
[CA i 2
B;:::;CB;D;0
s(f i)m
]
This equation is the starting pointfor all subspacealgorithms. An outline canbe given as
follows: 1. RegressY + t;f ontoZ t;p andU + t;f
inorder toobtainestimates ^ z and ^ u . 2. ^ z
willtypicallybeoffullrank,whereasthelimit
z isofrankn. ThusaSVDof ^ W + f ^ z ^ W p = ^ U ^ ^ V 0 = ^ U n ^ n ^ V 0 n + ^
Risperformedleadingtoanestimate ^ O f =( ^ W + f ) 1 ^ U n ; ^ K p = ^ n ^ V 0 n ( ^ W p ) 1 . Here ^ W + f and ^ W p
are weightingmatrices, which haveto be chosen by theuser. Wewill
assume throughout that they are nonsingular (a.s.) and converge to some limit W + f and W + p respectively forT !1. ^ n
contains then dominantsingular valuesand ^ U n ; ^ V n the
correspondingleftorrightrespectivelysingularvectors.
3. Eithertheestimates ^ O f and ^ u ortheestimate ^ K p
isusedtoestimatethesystem.
In the last step, one hastwo choices: MOESPand one variant of N4SIDuse the structure in O
f
and
u
in order toestimate thesystemmatrices, whereasCCAandanothervariantof N4SIDuse
^
K
p
toestimatethestateasx^
t = ^ K p Z t;p
andthenusethesystemequations(1)in ordertoobtain
estimatesofthesystem.
Nowtherearetwopossiblewaystoinclude thepreprocessinginthisframework:
Regressy t andu t ontop t
andlettheresidualsbedenotedasy~
t andu~
t
. This preprocessed
datathencanbeusedinthesubspacemethodin orderto obtainanestimateofthesystem
(A;B;C ;D;K). Use y t and z t =[u 0 t ;p 0 t ] 0
asthe datafor thesubspacemethod to obtainan estimateof the
system(A;B;C ;D;K)ofthetransferfunction relatingu
t toy
t .
The aim of the next section is to show that these two procedures deliver essentially identical
Note, thatin thepresentsetuptheregressionin therststeptoobtaintheestimates ^ z ; ^ u will
havetotakeintoaccountthesingularityofthematrixoftheregressorsduetothemulticollinearity
introducedbytheprocessp
t . ForZ
t;p
note, thatin theactualcalculationonlythevectorsZ ;
t;p
areused fortheestimate ^ z . Here Z ; t;p
denotestheresidualsof Z
t;p regressedonto U + t;f . Using thenotationha t ;b t i=T 1 P T f t=p+1 a t b 0 t thisisequalto Z t;p hZ t;p ;U + t;f ihU + t;f ;U + t;f i y U + t;f Here y
denotes the Moore Penrose pseudo inverse. Due to the nature of the variables p
t , the
residuals in all coordinates corresponding to p
t
are zero. Note that this corresponds to using
the preprocessed data y~
t ;u~
t
up to the dierence of at most f +p terms in the calculation of
theregressionsforthepreprocessing. Thisdierenceisasymptoticallynegligibleunder theusual
assumptionsonf andptobeo((logT) a
)forsomea<1. Thereforethecolumnsoftheestimate
^ z ,whichcorrepondtoy t oru t
respectively,areequivalenttotheonesobtainedfrompreprocessed
data. Theremainingcolumnscorrespondtop
t
andtheseareessentiallyzero(i.e. zerouptoeects
duetoinitialvaluesintheregressions). Alsonotethatthecalculationofthepseudoinversecould
becircumventedby usingonlythose coordinatesof U +
t;f
,whichare linearilyindependent. These
can be found easily, since the dependence structure of p
t
is known. In fact, if the regression
is performed accordingto theformula given above,all coordinates corresponding to p
t+j ;j >0
simplycanbeomittedwithoutchangingtheresult.
IntheMOESPtypeofproceduresthisshowsthat theestimates ^
O
f
ofthetwoalternativeswill
beessentiallyidentical. Thisshowstheequivalenceoftheestimates ^
Aand ^
C. Fortheestimation
ofB andD weconsideraspecialalgorithmsimilarto theMOESPtypeofprocedures: LetU
f ( ^ O f )
denotetheestimateofU
f
,whereelementsofO
f
arereplacedbytheirrespectiveestimates. Note
thatU f ( ^ O f
)islinearinBandD. ThentheestimateofBandDcanbeobtainedfromminimising
k ^ O ? f hY + t;f ;U + t;f ihU + t;f ;U + t;f i y U f ( ^ O f ) k Fr Here ^ O ? f ^ O f =0and ^ O ? f
isoffullrowrank. Fromthisitfollowsfromtheblockmatrixinversion,
that the columns correspondingto u
t
are identical to therespectivecolumns obtainedfrom the
preprocesseddata y~
t ;u~
t
. Fromthis itfollowsthattheestimation ofthepartofB andD, which
corresponds to u
t
is estimated essentiallyidentical asif the preprocesseddata would havebeen
used. This shows, that in this case theestimates of the coordinatesof thesystem (A;B;C ;D)
corresponding to u
t
is essentiallythesameasobtainedfrom thepreprocessed data. Concerning
thepartdue top
t
howeverthesameresultdoesnotseemtohold,comparetheexamplegivenin
section5. Thisisamatteroffurtherresearch.
FortheLarimoretypeofapproachwehaveseenthattheestimate ^
z
isessentiallyidenticalto
theestimateobtainedfromthepreprocesseddatainthecoordinatescorrespondingtotheoutputor
theinputduetou
t
,theothercomponentsbeingzero. Thereforethesameholdstruefor ^
K
p
,noting
that fortheusual weightingmatrices ^ W p ; ^ W p
alsoonlythedata Z ; t;p and Y +; t;f is used. Here Y +; t;f isdenedanalogouslyto Z ; t;f
. Thesevectorsareessentiallyidentical tothecorresponding
vectorsformed from the preprocesseddata. As aconsequencethe stateestimate x^
t = ^ K p Z t;p is
almostidenticaltotheestimateobtainedfromthepreprocesseddataexceptfordierencesinthe
spacespannedbythecomponentsofp
t
. Notethatthenextstepinthisclassofalgorithmsistouse
thestate in aregression,wherethese deterministic componentsare included. This againshows,
thattheestimateofthesystemmatricescorrespondingtou
t
coincideuptoinitialconditionswith
theestimateobtainedfromthepreprocesseddata. Note,that thenonidentiability inthis setup
iscircumvented,asthedeterministiccomponentsinx^
t+1 andx^
t
convergetosomevectors x and x respectively. Letting u
denote thecoeÆcients ofthe deterministiccomponentsof z
t
then it
followsfromthestateequationthat
x =A x +B u
t
any nonidentiability problems by construction and the estimates for B
p
, the columns due to
the deterministic components, are consistent and asymptotically normal. After the estimation
thenonidentiabilityofthemodelstructurecould beused toobtainadierentnormalisationas
B
p
=0usingthestructureofthenonidentiability. Inbothcasestheasymptoticdistributionwill
bethesame,asifthecorrespondingestimateshavebeencalculatedusing thepreprocesseddata,
asisstraightforwardtosee. Finallywehaveobtainedthefollowingtheorem:
Theorem1 Lety
t
be generatedbya systemofthe form(1), where the noise isi.i.d. withnite
momentsup tofourth order. Let the input z
t
consist of twodierent components: z
t =[u 0 t ;p 0 t ] 0 , where u t
is aquasistationary sequence, which is persistently exciting of order f+p. Further p
t
iseither a constant term,a trend component ora periodic component with known frequency and
shift. Then theestimatesofG
i
(q),the transferfunctionrelatingu
t toy
t
obtainedbytheinclusion
of p
t
in the MOESP type of subspace methods and the estimates obtained by using the residuals
from a regression of y t andu t onto p t
as the data input for the subspace method without p
t are
asymptotically equivalent, i.e. the asymptoticdistributionof the estimationerrorare identical.
Iftheinputu
t
isanARMAprocessgeneratedbyi.i.d. whitenoisehavingauniformelybounded
spectrumandp=p(T)=o((logT) a );p(T) dlogT 2logj 0 j
,thenthesameholdsfortheLarimoretype
ofalgorithms. Inthiscasealsotheestimatesofthetransferfunctionrelatingthenoisetotheoutput
willbeasymptoticallyequivalent. Furthermorealsotheestimates oftheeectsofthe deterministic
termsareessentially the same.
Itshould benoted,that theassumptionsof thetheorem arebynomeans necessary. Infact,
much weakerconditionssuÆce. Fora discussionon the necessaryconditionsin the MOESPtype
ofalgorithms see(Bauer andJansson, 2000),forassumptions inamartingaleframeworkfor the
Larimoretypeofproceduressee(Baueretal.,1999).
Thesignicanceof thetheorem is that itgivestheusertwodierent but equivalentwaysto
deal with nonzero means, driftsand periodic componentsincluded in manytime series. Either
one canpreprocess the data and estimate and test for the inclusion of these terms beforehand
andthenusethepreprocesseddatainorder toobtainestimatesofthedynamicalsystems,which
arethemaingoalin manyapplications. Thealternativewayisto usethedeterministictermsas
additionalinputsin thesubspacemethod. Bothprocedureshaveadvantagesand disadvantages.
The big advantage for the preprocessing approach lies in the fact, that the testing procedures,
whetherthedetmerinisticcomponentsarecontainedinthedataathand,arewellestablishedand
implementedin many programs. This facilitatestheanalysis fortheuser. Ontheotherhand in
thiscasetheinclusionofthepreprocessingstagesinthecalculationoftheasymptoticvarianceof
thesubspaceestimatesispossible,although cumbersome.
Theinclusionofthetermsintheestimationprocedureiscomputationallysimple. Theoriginal
procedurescanbeused,iftheregressionsaredoneinarobustwayusingthepseudoinversesinthe
caseofmulticollinearities. Thedirectembeddingofthepreprocessingmakesthecalculationofthe
varianceeasier,althoughatthepresenttimenoprogramstocalculatetheasymptoticvarianceof
the subspaceprocedure exists, that could be used in an industrial context. Forthe MOESPtype
ofproceduresthederivation of theasymptoticdistribution oftheestimates ofthecolumns ofB
and D corresponding to p
t
seemsto becomplicated. HoweverfortheCCA typeofestimates the
derivation of these distributions is straightforward. Tests for more realistic hypotheses can be
obtainedfrom theasymptotic distribution,e.g. whether ornot thedeterministiccomponentsin
theoutputcanbeexplainedonlythroughthecomponentsoccuringintheinput.
5 Numerical Exmaple
Inthissectionwewillpresentsomesimpleexamples,whichillustratethetheorygivenin thelast
section. TherstexampleisaonestateSISOsystem,whichisdescribedbythefollowingmatrices:
MOESP nopre. 0.9183 1.1498 0.0599 pre. 0.4983 0.9996 0.0013 incl. 0.4985 1.0001 0.0013 corr. 0.9983 0.9980 0.9978 CCA nopre. 0.8886 1.2125 0.0710 pre. 0.4999 0.9994 0.0011 incl. 0.5002 0.9998 0.0011 corr. 0.9970 0.9974 0.9974
Table1: ThistableshowsthemeanvaluesoftheestimatesofA;B andD estimatedwithsample
sizeT =1000usingMOESPandCCAontheoriginaldata(nopre.),onthepreprocesseddata(pre.)
orincludingaconstantterm(incl.) inthecalculations. Alsothecorrelation (corr.) betweenthe
estimates(pre.) and(incl.) isgiven.
TheinputisGaussianwhitenoisewithmean1andvariance1. Thenoise"
t
ischosento bezero
meanGaussianwhitenoiseofvariance1. Theoutputy
t isequalto y t =(D+C(q A) 1 B)u t +(1+C(q A) 1 K)" t +10
Here q denotes the forward shift operator. 1000 data sets of dimension T = 1000 have been
generatedandforeachdatasetthesystemisestimated
using nopreprocessingandnoinclusionofaconstant,
meancorrecteddataand
using theoriginaldatawiththeadditionofaconstantinputvariable.
f =p=5hasbeenchosenfortheMOESPtypeofprocedureandf =5;p=15fortheLarimoretype
ofprocedure,inallthetrials. Aftertheestimationthesystemshavebeentransformedtoechelon
canonicalform. InbothcasestheCCAweightingshavebeenused. Themeanofthecorresponding
estimatesfortheMOESPtypeofprocedures canbeseeninTable1.
Itcanbeseenclearly,thattheestimateswithouttakingthenonzeroconstantsintoaccountis
biased. It alsocanbeseen,thattheothertwoapproachesgivecomparablemeansandarehighly
correlated,asisexpectedfromthelastsection.
Comparingthemean ofy
t
accordingto theestimated model tothe actuallyestimated mean
on the other hand shows a dierent picture: For the CCA procedure a correlation of 0:9989 is
obtained,i.e. almostperfectagreement,whilefortheMOESPprocedurethecorrelationisestimated
as0:0948showingalmost no lineardependence. This reaÆrmsthe complicationsto be awaited
forthederivation oftheasymptoticdistribution ofthe estimatesofthecomponentsofB and D
correspondingtop
t .
Asasecondexamplethesamesystemisused,buttheinputischangedbyaddingarandomly
weighteddeterminsticcomponentandalsotheoutputiscontaminatedbyasimilarterm:
u t = v t +x t x t = 2 6 6 6 6 6 6 4 t sin(0:5t) sin(0:22t) 1 sin(0:5(t 1)) sin(0:22(t 1)) 3 7 7 7 7 7 7 5
0
20
40
60
80
100
120
140
160
180
200
−5
0
5
10
15
OUTPUT #1
0
20
40
60
80
100
120
140
160
180
200
−1
0
1
2
3
4
INPUT #1
Figure1: Observations1to200ofatypicaltrajecory.
Hereeachcomponentof isdrawnfroma[0;1]uniformdistribution. Also
y t =G(q)u t +H(q)" t +x t where G(q) = D+C(qI A) 1 B H(q) = I+C(qI A) 1 K
Thecomponentsof areuniformely[0;1]distributed. Finallyv
t and"
t
areGaussian zeromean
and unit variance random processes. 1000 time series with these characteristicsof sample size
T =1000areconstructed. Partsofatypicalinput-outputtrajectorypaircanbeseeninFigure1.
Theestimationshowsthesamepictureastherstexample:Figure2showsthestandarddeviations
ofthe estimatesofthe transferfunction D+C(qI A) 1
B at50 pointsin theanguar frequency
range[;]. Herebyf =p=5wasused andtheMOESPtypeof procedures. Theplotshowsthe
standarddeviationoftheMOESPestimatesusingthepreprocesseddataandthestandarddeviation
of the dierence between the estimates obtained from the preprocessed data and the estimates
obtained using the original data and including the additional variables, respectively. It canbe
seen clearly, that the dierence in between the twoapproaches is of lowermagnitude than the
estimation errors themselves. For the CCA case the results are almost identical, therefore the
presentationoftheresultsisomitted.
6 Conclusions
In this paper we did derive the asymptotic properties of subspace procedures, when the data,
whichisusedfortheprocedure,ispreprocessedbyremovingtrendsandperiodiccomponents. It
has been shown, that the preprocessingcan alternatively beseen asthe inclusion of additional
inputterms. Thismakesitpossibletocalculatetheasymptoticvarianceoftheestimatesobtained
withtheprocedurefromthestandardproceduresforthesubspaceestimates. Alsotestprocedures
−1
0
1
10
−5
10
−4
10
−3
10
−2
10
−1
10
0
Ang. frequ. /
π
Standard Deviation
Comparison of Standard Deviations
MOESP
Diff.
Figure2: Theplotsshowthestandarddeviationofthetransferfunctionestimatesobtainedfrom
thepreprocesseddataandf =p=5(|) andthestandarddeviationofthedierenceofthetwo
dierentprocedures(-+-).
References
Bauer,D.(1998).SomeAsymptoticTheoryfortheEstimationofLinearSystemsUsingMaximum
LikelihoodMethodsorSubspaceAlgorithms.PhDthesis.TUWien.
Bauer, D. and M.Jansson (2000).Analysis of theasymptoticpropertiesof theMOESP typeof
subspacealgorithms.Automatica36(4),497{509.
Bauer, D., M. Deistler and W. Scherrer (1999). Consistency and asymptotic normality of some
subspacealgorithmsforsystemswithoutobservedinputs.Automatica35,1243{1254.
Larimore,W. E. (1983).System identication, reduced order ltersand modeling viacanonical
variateanalysis.In: Proc.1983Amer.ControlConference2.(H.S.RaoandP.Dorato,Eds.).
Piscataway,NJ.pp.445{451. IEEEServiceCenter.
Van Overschee, P. and B. DeMoor(1994). N4sid: Subspace algorithms for the identication of
combineddeterministic-stochasticsystems.Automatica30, 75{93.
Verhaegen,M. (1994).Identication of thedeterministic partofmimo statespace modelsgiven