The question makes little or no sense except in projected coordinates, because a "bounding box" and even a constant bearing have ambiguous meanings on a spheroid. In the projected plane, all you have to do is rotate the points to make the bearing horizontal or vertical. Finding the bounding box is then done the usual way (by obtaining the extreme values of the coordinates). Rotating this figure back finishes the job.
To maintain high precision, rotate the points around some central (or nearby) location.
If we were to imagine a GIS that supported the basic operations of (a) rotating figures and (b) finding bounding boxes--which most of them do--its solution would read like this:
BoundingBox = Rotate(Extent(Rotate(points, bearing)), -bearing)
where
Rotate(p, a)
rotates a feature p
by an amount a
and
Extent(p)
returns the (rectangular) extent of a feature.
Here is an example showing the rotated and original situations:

R
code (which is readily ported to Python or any other platform supporting basic matrix operations) follows. Most of it just generates sample data and plots the results.
#
# Sample data (looking like those of the question,
# with typical ranges of projected coordinates).
#
xy <- cbind(1:4 + 200000, c(0,5,0,4) + 4000000)
bearing <- 45 * pi/180 # Radians east of north
#
# Helper functions.
#
rotate <- function(xy, a, origin=c(0,0)) {
c <- cos(a); s <- sin(a)
return(t(matrix(c(c, -s, s, c), 2) %*% (t(xy) - origin) + origin))
}
extent <- function(xy) {
e <- apply(xy, 2, range)
return(matrix(c(e[1,1], e[2,1], e[2,1], e[1,1],
e[1,2], e[1,2], e[2,2], e[2,2]), ncol=2, byrow=FALSE))
}
#
# Compute the oriented bounding box.
#
center <- apply(xy, 2, mean)
bb <- rotate(extent(rotate(xy, bearing, center)), -bearing, center)
#
# Display the points and their oriented bounding box.
#
plot(rbind(bb, xy), type="n", asp=1, xlab="X", ylab="Y", main="Solution")
polygon(bb, col="#f0f0f0")
points(xy, pch=19, col="Red")