Problem simulating AR(2) process

Question

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.

Answer 1

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...