## 2 Brief Review of MIMO-NARX Modeling  
### MIMO-NARX Model for MDOF System  
In structural dynamics, a nonlinear MDOF dynamic system can be expressed as:  
\[\mathbf{M}\ddot{\mathbf{x}}(t)+\mathbf{C}\dot{\mathbf{x}}(t)+\mathbf{f}_{r}(\mathbf{x}(t))=\mathbf{f}(t) \tag{1}\]  
where \(\mathbf{M}\) and \(\mathbf{C}\) are the mass and damping matrix, respectively; \(t\) denotes time; \(\ddot{\mathbf{x}}(t)\), \(\dot{\mathbf{x}}(t)\) and \(\mathbf{x}(t)\) are response vectors of acceleration, velocity, and displacement, respectively; \(\mathbf{f}_{r}\) is the vector of nonlinear hysteresis forces; \(\mathbf{f}(t)\) is the external excitation force, which is often expressed in earthquake engineering as \(\mathbf{f}(t)=-\mathbf{M}\ddot{\mathbf{x}}_{g}(t)\) with an all-ones vector \(\mathbf{\iota}\) and the seismic acceleration \(\ddot{\mathbf{x}}_{g}(t)\). The nonlinear dynamic behavior of the MDOF system in Equation (1) can be captured in a NARX model by considering several steps of the current and past responses of input and output. A MIMO-NARX technique is adopted in this study to capture coupled nonlinear dynamics since a SISO-NARX model is unable to represent the response coupling between multiple outputs for general nonlinearity [13]. The current output \(\mathbf{x}(t)\) of nonlinear MDOF dynamic system depend on its past output values \([\mathbf{x}(t-\Delta t),\ldots,\mathbf{x}(t-n_{o}\Delta t)]\) with maximum time lag \(n_{o}\) and time interval \(\Delta t\), and current and past input values \([\mathbf{f}(t),\mathbf{f}(t-\Delta t),\ldots,\mathbf{f}(t-n_{i}\Delta t)]\) with maximum time lag \(n_{i}\) and can be expressed as:  
\[\mathbf{x}(t)=\mathbf{G}(\mathbf{z}(t))+\mathbf{\epsilon}(t) \tag{2}\]  
where \(\mathbf{z}(t)=\left[\mathbf{x}^{\rm T}(t-\Delta t),\ldots,\mathbf{x}^{\rm T}(t-n_{o} \Delta t),\mathbf{f}^{\rm T}(t),\ldots,\mathbf{f}^{\rm T}(t-n_{i}\Delta t)\right]\) is the regression vector of current and past values of input and output; \(\mathbf{G}(\cdot)\) represents the MIMO-NARX model to be identified based on the input and output; and \(\mathbf{\epsilon}(t)\) is vector-valued Gaussian error process [20]. This study adopts the linear-in-the-parameter form for the MIMO-NARX model which can be expressed as:  
\[\mathbf{G}(\mathbf{z}(t))=\mathbf{\Theta}^{\rm T}\mathbf{g}(\mathbf{z}(t)) \tag{3}\]  
where \(\mathbf{g}(\cdot)=\left[\mathbf{g}_{1}^{\rm T}(\cdot);\mathbf{g}_{2}^{\rm T}(\cdot), \ldots,\mathbf{g}_{n}^{\rm T}(\cdot)\right]^{\rm T}\) contains the NARX terms for the \(n\) DOFs system; and \(\mathbf{\Theta}=\text{diag}[\mathbf{\Theta}_{1},\ldots,\mathbf{\Theta}_{n}]\) contains the NARX coefficients of the \(n\) NARX models. For each degree of freedom, the output of a MIMO-NARX model can be formulated as:  
\[x_{j}(t)=\mathbf{\Theta}_{j}^{\rm T}\mathbf{g}_{j}(\mathbf{z}(t))+\mathbf{\epsilon}_{j}(t) \tag{4}\]  
where \(j\) represents the \(j\)th degree of freedom; \(\mathbf{g}_{j}(\mathbf{z}(t_{i}))\) is the most relevant NARX terms selected from all potential NARX terms, that is, predetermined basis functions. While only the \(j\)th degree of freedom is considered in the SISO-NARX strategy, the MIMO-NARX strategy includes the input and output histories from all degrees of freedom. Identification of the MIMO-NARX model requires the selection of NARX terms and calibration of their coefficient. An efficient approach [18] has been proposed for identifying the NARX model for SISO systems by implementing the least-angle regression (LARs) algorithm for structure determination and the ordinary least square (OLS) method for coefficient calibration. The LARs algorithm [21] is employed to select candidate NARX model terms by calculating the correlation between each potential model term and the system output [13].