In the Time for action – filtering a detrended signal section we created a simple filter for detrended data. Now let’s use a more restrictive filter that will leave us only with the main frequency component. We will fit a sinusoidal pattern to it and plot our results. This model has four parameters—amplitude, frequency, phase, and vertical offset. Perform the following steps to fit to a sine:
- Define a
residualsfunction based on a sine wave model.def residuals(p, y, x): A,k,theta,b = p err = y-A * np.sin(2* np.pi* k * x + theta) + b return err - Transform the filtered signal back to the original domain.
- Guess the values of the parameters for which we are trying to estimate a transformation from the time domain into the frequency domain.
N = len(qqq) f = np.linspace(-N/2, N/2, N) p0 = [filtered.max(), f[amps.argmax()]/(2*N), 0, 0] print “P0”, p0
The initial values would be shown as follows:
P0 [2.6679532410065212, 0.00099598469163686377, 0, 0] - Call the
leastsqfunction.plsq = optimize.leastsq(residuals, p0, args=(filtered, dates)) p = plsq[0] print “P”, p
The following are the final parameter values:
P [ 2.67678014e+00 2.73033206e-03 -8.00007036e+03 -5.01260321e-03] - Finish the first subplot with detrended data, filtered data, and fit of the filtered data. Use a date format for the horizontal axis and add a legend.
plt.plot(dates, y, ‘o’, label=”detrended”) plt.plot(dates, filtered, label=”filtered”) plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], ‘^’, label=”fit”) fig.autofmt_xdate() plt.legend(prop={‘size’:’x-large’}) - Add a second subplot with a legend of the main component of the frequency spectrum.
ax2 = fig.add_subplot(212) plt.plot(f, amps, label=”transformed”)
The following shows the resulting charts:
We detrended 1 year of price data for QQQ. This signal was then filtered until only the main component of the frequency spectrum was left over. We fitted a sine to the filtered signal using the
scipy.optimize module (see optfit.py).
from matplotlib.finance import quotes_historical_yahoo
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
from scipy import signal
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
from scipy import optimize
start = (2010, 7, 25)
end = (2011, 7, 25)
quotes = quotes_historical_yahoo(“QQQ”, start, end)
quotes = np.array(quotes)
dates = quotes.T[0]
qqq = quotes.T[4]
y = signal.detrend(qqq)
alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter(“%b %Y”)
fig = plt.figure()
fig.subplots_adjust(hspace=.3)
ax = fig.add_subplot(211)
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(month_formatter)
ax.tick_params(axis=’both’, which=’major’, labelsize=’x-large’)
amps = np.abs(fftpack.fftshift(fftpack.rfft(y)))
amps[amps < amps.max()] = 0
def residuals(p, y, x):
A,k,theta,b = p
err = y-A * np.sin(2* np.pi* k * x + theta) + b
return err
filtered = -fftpack.irfft(fftpack.ifftshift(amps))
N = len(qqq)
f = np.linspace(-N/2, N/2, N)
p0 = [filtered.max(), f[amps.argmax()]/(2*N), 0, 0]
print “P0”, p0
plsq = optimize.leastsq(residuals, p0, args=(filtered, dates))
p = plsq[0]
print “P”, p
plt.plot(dates, y, ‘o’, label=”detrended”)
plt.plot(dates, filtered, label=”filtered”)
plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], ‘^’, label=”fit”)
fig.autofmt_xdate()
plt.legend(prop={‘size’:’x-large’})
ax2 = fig.add_subplot(212)
ax2.tick_params(axis=’both’, which=’major’, labelsize=’x-large’)
plt.plot(f, amps, label=”transformed”)
plt.legend(prop={‘size’:’x-large’})
plt.show()