Bayesian point estimation is commonly used for system identification owing to its good properties for small sample sizes. Although this type of estimator is usually non-parametric, Bayes estimates can also be obtained for rational parametric models, which is often of interest. However, as in maximum-likelihood methods, the Bayes estimate is typically computed via local numerical optimization that requires good initialization and cannot guarantee global optimality. In this contribution, we propose a computationally tractable method that computes the Bayesian parameter estimates with posterior certification of global optimality via sum-of-squares polynomials and sparse semidefinite relaxations. It is shown that the method is applicable to certain discrete-time linear models, which takes advantage of the rational structure of these models and the sparsity in the Bayesian parameter estimation problem. The method is illustrated on a simulation model of a resonant system that is difficult to handle when the sample size is small.