This paper presents a method to extend the system identification method for multivariate time series (e.g., neural and behavioral records) by applying preferred subspace identification (PSID) for filtering and smoothing. While the existing PSID only uses the historical main time series data to predict the auxiliary time series, this paper enables better estimation by incorporating the concurrent data (filtering) or all the data (smoothing). First, we present a method to uniquely identify the model with the optimal Kalman update step by exploiting the presence of auxiliary signals, and for this purpose, we develop PSID filtering with an additional reduced rank regression step. Second, inspired by the two-filter Kalman smoother formulation, we develop a forward-backward PSID smoothing algorithm to reapply PSID filtering backwardly on the residual of the filtered auxiliary signals. We verify through simulation data that the proposed method recovers the true model parameters for filtering and achieves optimal filtering and smoothing decoding performance that matches the ideal performance of the actual underlying model. In conclusion, this study provides a principled framework for optimal linear filtering and smoothing in two-signal settings, significantly expanding the toolkit for dynamic interaction analysis of multivariate time series.