Result
First we state the result, then we will go on to derive and explain it.

= 


where:
 r_{00} = cos(θ), r_{01} = sin(θ), r_{10} = sin(θ), r_{11} = cos(θ)
 θ = angle we are rotating around (in appropriate units for the trig functions we are using)
 x_{in},y_{in} = the coordinates of the input point
 x_{out},y_{out} = the coordinates of the output point (the result)
 x,y = the coordinates of the point that we are rotating around.
Or multiplying out the matrix and vector terms to give in ordinary equations:
x_{out} = r_{00}* x_{in} + r_{01}* y_{in} + x  r_{00}*x  r_{01}*y
y_{out} = r_{10}* x_{in} + r_{11}* y_{in} + y  r_{10}*x  r_{11}*y
Explanation
In order to calculate the rotation about any arbitrary point we need to calculate its new rotation and translation. In other words rotation about a point is an 'proper' isometry transformation' which means that it has a linear and a rotational component.
Assume we have a matrix [R0] which defines a rotation about the origin:
We now want to apply this same rotation but about an arbitrary point P:
As we can see its orientation is the same as if it had been rotated about the origin, but it has been translated to a different point on space by the rotation.
In order to prove this and to calculate the amount of linear translation we need to replace:
 translate about arbitrary point P (Px,Py).
With the following 3 simpler transforms which, when done in order, are equivalent:
 translate the arbitrary point to the origin (subtract P which is translate by Px,Py)
 rotate about the origin (can use 2×2 matrix R0)
 then translate back. (add P which is translate by +Px,+Py)
So if we are using the global frameofreference (as explained here)
then,
[resulting transform] = [third transform] * [second transform] * [first transform]
[resulting transform] = [+Px,+Py] * [rotation] * [Px,Py]
Note for matrix algebra, the order of operations is important, so these translations do not cancel out.
So matrix representing rotation about a given point is:
[R] = [T]^{1} * [R0] * [T]
where:
[T]^{1} = inverse transform = translation of point to origin
1  0  x 
0  1  y 
0  0  1 
[R0] = rotation about origin (if this is not clear see this discussion)
r_{00} = cos(θ)  r_{01} = sin(θ)  0 
r_{10} = sin(θ)  r_{11} = cos(θ)  0 
0  0  1 
[T] = translation of origin to point
1  0  x 
0  1  y 
0  0  1 
when these matrices are multiplied this will give the following result for rotation about x,y:

* 

* 

multiply out second pair

* 

multiply remaining pair:
r_{00}  r_{01}  x  r_{00}*x  r_{01}*y 
r_{10}  r_{11}  y  r_{10}*x  r_{11}*y 
0  0  1 
So the rotational components are the same but the rotation moves the position of the centre.
Further Reading
You may be interested in other means to represent orientation and rotational quantities such as:
Or you may be interested in how these quantities are used to simulate physical objects: