3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – vinux Feb 13 '12 at 18:31
  • $\begingroup$ 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)? $\endgroup$ – Xi'an Feb 13 '12 at 19:42
  • 1
    $\begingroup$ @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]. $\endgroup$ – Néstor Feb 13 '12 at 19:54
  • $\begingroup$ Then I do not see the reason why it fails... $\endgroup$ – Xi'an Feb 13 '12 at 20:55
  • $\begingroup$ Nice to see you solved the problem. $\endgroup$ – Xi'an Feb 13 '12 at 21:29
3
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.