While subspace identification methods (SIMs) are appealing due to their simple parameterization for MIMO systems and robust numerical realizations, a comprehensive statistical analysis of SIMs remains an open problem, especially in the non-asymptotic regime. In this work, we provide a finite sample analysis for a class of SIMs, which reveals that the convergence rates for estimating Markov parameters and system matrices are O(1 / SN), in line with classical asymptotic results. Based on the observation that the model format in classical SIMs is non-causal because of a projection step, we choose a parsimonious SIM that bypasses the projection step and strictly enforces a causal model to facilitate the analysis, where a bank of ARX models are estimated in parallel. Leveraging recent results from a finite sample analysis of an individual ARX model, we obtain a union error bound for an array of ARX models and proceed to derive error bounds for system matrices using robustness results for the singular value decomposition.