Assuming weak stationarity of AR process:
In[1]:= Variance[ARProcess[{a1, a2}, var][t]]
Out[1]= -(((-1 + a2) var)/((1 + a2) (1 - a1^2 - 2 a2 + a2^2)))
gives the answer for p=2 and for var denoting the variance of white noise driving the process. In Mathematica ARProcess is assumed to be driven by a Gaussian white noise, but this formula is more general - assuming zero first moment and finite second moment of the noise. We also assume that the noise variable is independent of past observations.
Clearly it becomes more complicated in higher dimensions:
In[2]:= Variance[ARProcess[{Array[a, {2, 2}]}, Array[v, {2, 2}]][t]]
Out[2]= {-((v[1, 1] - a[1, 2] a[2, 1] v[1, 1] -
a[1, 1] a[2, 2] v[1, 1] - a[2, 2]^2 v[1, 1] -
a[1, 2] a[2, 1] a[2, 2]^2 v[1, 1] + a[1, 1] a[2, 2]^3 v[1, 1] +
a[1, 1] a[1, 2] v[1, 2] + a[1, 2]^2 a[2, 1] a[2, 2] v[1, 2] -
a[1, 1] a[1, 2] a[2, 2]^2 v[1, 2] + a[1, 1] a[1, 2] v[2, 1] +
a[1, 2]^2 a[2, 1] a[2, 2] v[2, 1] -
a[1, 1] a[1, 2] a[2, 2]^2 v[2, 1] + a[1, 2]^2 v[2, 2] -
a[1, 2]^3 a[2, 1] v[2, 2] +
a[1, 1] a[1, 2]^2 a[2, 2] v[2, 2])/(-1 + a[1, 1]^2 +
a[1, 2] a[2, 1] + a[1, 1]^2 a[1, 2] a[2, 1] +
a[1, 2]^2 a[2, 1]^2 - a[1, 2]^3 a[2, 1]^3 + a[1, 1] a[2, 2] -
a[1, 1]^3 a[2, 2] + 3 a[1, 1] a[1, 2]^2 a[2, 1]^2 a[2, 2] +
a[2, 2]^2 - a[1, 1]^2 a[2, 2]^2 + a[1, 2] a[2, 1] a[2, 2]^2 -
3 a[1, 1]^2 a[1, 2] a[2, 1] a[2, 2]^2 - a[1, 1] a[2, 2]^3 +
a[1, 1]^3 a[2, 2]^3)), -((-a[2, 1]^2 v[1, 1] +
a[1, 2] a[2, 1]^3 v[1, 1] - a[1, 1] a[2, 1]^2 a[2, 2] v[1, 1] -
a[1, 1] a[1, 2] a[2, 1]^2 v[1, 2] - a[2, 1] a[2, 2] v[1, 2] +
a[1, 1]^2 a[2, 1] a[2, 2] v[1, 2] -
a[1, 1] a[1, 2] a[2, 1]^2 v[2, 1] - a[2, 1] a[2, 2] v[2, 1] +
a[1, 1]^2 a[2, 1] a[2, 2] v[2, 1] - v[2, 2] +
a[1, 1]^2 v[2, 2] + a[1, 2] a[2, 1] v[2, 2] +
a[1, 1]^2 a[1, 2] a[2, 1] v[2, 2] + a[1, 1] a[2, 2] v[2, 2] -
a[1, 1]^3 a[2, 2] v[2, 2])/(1 - a[1, 1]^2 - a[1, 2] a[2, 1] -
a[1, 1]^2 a[1, 2] a[2, 1] - a[1, 2]^2 a[2, 1]^2 +
a[1, 2]^3 a[2, 1]^3 - a[1, 1] a[2, 2] + a[1, 1]^3 a[2, 2] -
3 a[1, 1] a[1, 2]^2 a[2, 1]^2 a[2, 2] - a[2, 2]^2 +
a[1, 1]^2 a[2, 2]^2 - a[1, 2] a[2, 1] a[2, 2]^2 +
3 a[1, 1]^2 a[1, 2] a[2, 1] a[2, 2]^2 + a[1, 1] a[2, 2]^3 -
a[1, 1]^3 a[2, 2]^3))}