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=W X=W X=a1*X+a2*X+W X=a1*X+a2*X+W ... 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.