Your explanation is wonderful!! I think this is exactly what I needed to know to work my setting which is based on this. Thank you so much!!!!!

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))}