Given massive data that may be time dependent and multi-dimensional, how to efficiently explore the underlying functional relationships across different dimensions and time lags? In this work, we propose a methodology to sequentially and adaptively model nonlinear multivariate time series data. Data at each time step and dimension is modeled as a nonlinear function of past values corrupted by noise, and the underlying nonlinear function is assumed to be approximately expandable in a spline basis. We cast the modeling of data as finding a good fit representation in the linear span of multi-dimensional spline basis, and use a variant of h-penalty regularization in order to reduce the dimensionality of representation. Using adaptive filtering techniques, we design our online algorithm to automatically tune the underlying parameters based on the minimization of the regularized sequential prediction error. We demonstrate the generality and flexibility of the proposed approach on both synthetic and real-world datasets. Moreover, we analytically investigate the performance of our algorithm by obtaining bounds of the prediction errors.