The intercept is already included in the original model: see the matrix of covariates X contains a columns of 1s.
In your code above, since the linear regression model only contains a intercept mu, you can compute the eta = np.exp(trace['mu']) and beta = 1/trace['s']