|
| 1 | +from matplotlib.finance import quotes_historical_yahoo |
| 2 | +import numpy as np |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +from scipy import fftpack |
| 5 | +from scipy import signal |
| 6 | +from matplotlib.dates import DateFormatter |
| 7 | +from matplotlib.dates import DayLocator |
| 8 | +from matplotlib.dates import MonthLocator |
| 9 | +from scipy import optimize |
| 10 | + |
| 11 | + |
| 12 | +start = (2010, 7, 25) |
| 13 | +end = (2011, 7, 25) |
| 14 | + |
| 15 | +quotes = quotes_historical_yahoo("QQQ", start, end) |
| 16 | +quotes = np.array(quotes) |
| 17 | + |
| 18 | +dates = quotes.T[0] |
| 19 | +qqq = quotes.T[4] |
| 20 | + |
| 21 | + |
| 22 | +y = signal.detrend(qqq) |
| 23 | + |
| 24 | + |
| 25 | +alldays = DayLocator() |
| 26 | +months = MonthLocator() |
| 27 | +month_formatter = DateFormatter("%b %Y") |
| 28 | + |
| 29 | +fig = plt.figure() |
| 30 | +fig.subplots_adjust(hspace=.3) |
| 31 | +ax = fig.add_subplot(211) |
| 32 | + |
| 33 | +ax.xaxis.set_minor_locator(alldays) |
| 34 | +ax.xaxis.set_major_locator(months) |
| 35 | +ax.xaxis.set_major_formatter(month_formatter) |
| 36 | +ax.tick_params(axis='both', which='major', labelsize='x-large') |
| 37 | + |
| 38 | +amps = np.abs(fftpack.fftshift(fftpack.rfft(y))) |
| 39 | +amps[amps < amps.max()] = 0 |
| 40 | + |
| 41 | +def residuals(p, y, x): |
| 42 | + A,k,theta,b = p |
| 43 | + err = y-A * np.sin(2* np.pi* k * x + theta) + b |
| 44 | + |
| 45 | + return err |
| 46 | + |
| 47 | +filtered = -fftpack.irfft(fftpack.ifftshift(amps)) |
| 48 | +N = len(qqq) |
| 49 | +f = np.linspace(-N/2, N/2, N) |
| 50 | +p0 = [filtered.max(), f[amps.argmax()]/(2*N), 0, 0] |
| 51 | +print "P0", p0 |
| 52 | + |
| 53 | +plsq = optimize.leastsq(residuals, p0, args=(filtered, dates)) |
| 54 | +p = plsq[0] |
| 55 | +print "P", p |
| 56 | +plt.plot(dates, y, 'o', label="detrended") |
| 57 | +plt.plot(dates, filtered, label="filtered") |
| 58 | +plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], '^', label="fit") |
| 59 | +fig.autofmt_xdate() |
| 60 | +plt.legend(prop={'size':'x-large'}) |
| 61 | + |
| 62 | +ax2 = fig.add_subplot(212) |
| 63 | +ax2.tick_params(axis='both', which='major', labelsize='x-large') |
| 64 | +plt.plot(f, amps, label="transformed") |
| 65 | + |
| 66 | +plt.legend(prop={'size':'x-large'}) |
| 67 | +plt.show() |
0 commit comments