# Problem simulating AR(2) process

OP EDIT: There where no problem with this. The problem was with the method I was using for obtaining the PACF. Apparently it doesn't work quite well in this case (I was using the scikits/tsa python package to obtain the PACF via the YW equations). Testing the coefficients in R worked like a charm.

I'm trying to simulate an AR(2) process, but I seem to be getting awful results.

The way I'm doing it is as follows: if I want to simulate 1000 points of an AR(2) process with coefficients, say, a1=0.1 and a2=0.5, I simulate a realization of 2000 points of a white noise process (in my case I simulated 2000 points drawn from a normal, ~N(0,1) distribution), where I'll use the first 1000 points as burn-in points. Suppose I store this realization in a vector W[t]. Then, I simulate the actual AR(2) process by iterating on a new vector, X[t], as follows:

X[0]=W[0]

X[1]=W[1]

X[2]=a1*X[1]+a2*X[0]+W[2]

X[3]=a1*X[2]+a2*X[1]+W[3]

...

X[i]=a1*X[i-1]+a2*X[i-2]+W[i]


Finally, I return the last 1000 values of the vector X[i]. The problem is that when I plot the actual Partial Autocorrelation Function (i.e. when I estimate the coefficients of the AR(2) process I generated), I get wrong coefficients for a1 (I get an acceptable value for coefficient a2, though). What am I missing?

Here's my Python function for the simulation:

from numpy.random import normal
from pylab import *

# INPUT:
# a:        Is the array with coefficients, e.g. a=array([a1,a2]).
# sigma:    The white noise (zero-mean normal in this case) standard deviation.
# n:        Number of points to generate.

def ARgenerator(a,sigma,n,burnin=0):

if(burnin==0):
burnin=100*len(a) # Burn-in elements!
w=normal(0,sigma,n+burnin)
AR=array([])
s=0.0
warning=0
for i in range(n+burnin):
if(i<len(a)):
AR=append(AR,w[i])
else:
s=0.0
for j in range(len(a)):
s=s+a[j]*AR[i-j-1]
AR=append(AR,s+w[i])
print 'Measured standard deviation: '+str(sqrt(var(w[burnin:])))
return AR[burnin:]


Edit: by an MCMC simulation (where I use the autocovariance matrix for an AR(2) process with a multi-variate gaussian likelihood), I get right the value of sigma (the std. deviation of the white noise process) and the a2 coefficient. However, the a1 value I obtain has nothing to do with the true one.

-
Logic wise there is no issue. I would suggest you to do the same by changing the random seed. If you get the same problem, this probably due to some program error. Without the program or data, I cannot really answer what is the problem. –  vinux Feb 13 '12 at 18:31
I do not know Python and the following is likely to be silly, but isn't there a problem with the (i-j-1) as the index of AR. Shouldn't it be (i-j)? –  Xi'an Feb 13 '12 at 19:42
@Xi'an: No. For an array of coefficients of length 2, len(a)=2. Therefore, the line "for j in range(len(a))" iterates j=0 and j=1. On the other hand, the iteration is done only if i>=2, so in the first iteration i=2. Therefore, s=s+a[j]*AR[i-j-1] for j=0 is s=s+a[0]*AR[1], and for j=1 is s=s+a[1]*AR[0]. –  Néstor Feb 13 '12 at 19:54
Then I do not see the reason why it fails... –  Xi'an Feb 13 '12 at 20:55
Nice to see you solved the problem. –  Xi'an Feb 13 '12 at 21:29

This seems right, so there may be a mistake in your actual simulation... Here is the R version:

    > eps=rnorm(10^3)
> x=rnorm(10^3)
> a=0.1;b=0.5
> x[2]=a*x[1]+eps[2]
> for (t in 3:10^3) x[t]=a*x[t-1]+b*x[t-2]+eps[t]
> arima(x,c(2,0,0),incl=FALSE)

Call:
arima(x = x, order = c(2, 0, 0), include.mean=FALSE)

Coefficients:
ar1     ar2
0.0859  0.5001
s.e.  0.0274  0.0274

sigma^2 estimated as 0.9547:  log likelihood = -1396.07,  aic = 2800.14


which clearly recovers the coefficients...

-
I added my Python code. I can't seem to get where my mistake is. –  Néstor Feb 13 '12 at 19:18