[Numpy-discussion] Recurrence relationships

[josef.pktd at gmail.com]

On Wed, May 6, 2009 at 10:21 AM,  <josef.pktd at gmail.com> wrote:
> On Wed, May 6, 2009 at 10:00 AM, Talbot, Gerry <Gerry.Talbot at amd.com> wrote:
>> Sorry, I guess I wasn't clear, I meant:
>>
>>        for n in xrange(1,N):
>>          y[n] = A*x[n] + B*y[n-1]
>>
>> So y[n-1] is the result from the previous loop iteration.
>>
>
> I was using scipy.signal for this but I have to look up what I did
> exactly. I think either signal.correlate or using signal.lti.
>

No, its signal.lfilter, below is a part of a script I used to simulate
and estimate an AR(1) process, which is similar to your example.

I haven't looked at it in a while but it might give you the general idea.

Josef

# Simulate AR(1)
#--------------
# ar * y = ma * eta
ar = [1, -0.8]
ma = [1.0]

# generate AR data
eta = 0.1 * np.random.randn(1000)
yar1 = signal.lfilter(ar, ma, eta)

etahat = signal.lfilter(ma, ar, y)
np.all(etahat == eta)

# find error for given filter on data
print 'AR(2)'
for rho in [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.79, 0.8, 0.81, 0.9]:
    etahatr = signal.lfilter(ma, [1, --rho], yar1)
    print rho,np.sum(etahatr*etahatr)
print 'AR(2)'
for rho2 in np.linspace(-0.4,0.4,9):
    etahatr = signal.lfilter(ma, [1, -0.8, -rho2], yar1)
    print rho2,np.sum(etahatr*etahatr)

def errfn(rho):
    etahatr = signal.lfilter(ma, [1, -rho], yar1)
    #print rho,np.sum(etahatr*etahatr)
    return etahatr

def errssfn(rho):
    etahatr = signal.lfilter(ma, [1, -rho], yar1)
    return np.sum(etahatr*etahatr)


resultls = optimize.leastsq(errfn,[0.5])
print 'LS ARMA(1,0)', resultls

resultfmin = optimize.fmin(errssfn, 0.5)
print 'fminLS ARMA(1,0)', resultfmin