Python Pandas Kaufman自适应移动平均(KAMA) - Pandas或Cython中的递归计算

Question

我正在尝试在Python Pandas或Cython中为考夫曼自适应移动平均(KAMA)创建一个函数(我已经在R&Rcpp中做到了这一点)。我有关于递归计算的问题，filt1

filt1[i] =  ( filt1[i-1] + SC[i]*(price[i]-filt1[i-1]) )

我期望KAMA系列应该有(i) NA's开始，n=10的长度(ii)开始KAMA，对于2010年1月19日，价格的原始平均值1142.393，在这种情况下，收盘价的平均值(iii)之后的KAMA值来自递归公式filt1i

所以：

KAMA
2010-01-04       NA
2010-01-05       NA
2010-01-06       NA
2010-01-07       NA
2010-01-08       NA
2010-01-11       NA
2010-01-12       NA  
2010-01-13       NA
2010-01-14       NA
2010-01-15       NA
2010-01-19 1142.393
2010-01-20 1142.367
2010-01-21 1142.244
2010-01-22 1140.212
2010-01-25 1138.683
2010-01-26 1136.517

以…结尾

2013-12-24 1791.114
2013-12-26 1802.816
2013-12-27 1814.759
2013-12-30 1822.844
2013-12-31 1830.523

我已经开始了，

#%%
# Include this line for NEW WINDOW(S) for figures
%pylab qt4 

# start with getting some data to test on
import datetime
import tradingWithPython as twp # main toolkit functions
import tradingWithPython.lib.yahooFinance as yf # yahoo finance module   
import tradingWithPython.lib.backtest as backtest
from tradingWithPython.lib.extra import ProgressBar # import progress bar

#The python module (talib) that I will be using to calculate the technical 
#indicators is a wrapper around the open source TA-Lib. 
import talib
import numpy as np

import pandas.io.data as web  

import pandas.stats.moments

import pandas as pd
#pd.set_option('html', False) # how to display data - DEFAULT is True
#pd.set_option('display.height', int(1e7))
pd.set_option('display.max_rows', int(1e7))
#pd.set_option('display.max_columns', int(1e7))
pd.set_option('display.width', 3000)
#%%


#%%
def KAMA(x, n=10, pow1=2, pow2=30):
  ''' kama indicator '''    
  ''' accepts pandas dataframe of prices '''


  d['absDiffx'] = abs(x - x.shift(1) )  

  d['ER.num'] = ( x - x.shift(n) )
  d['ER.den'] = pandas.stats.moments.rolling_sum(d['absDiffx'],n)
  d['ER'] = d['ER.num'] / d['ER.den']

  d['SC'] = ( d['ER']*(2.0/(pow1+1)-2.0/(pow2+1.0))+2/(pow2+1.0) ) ** 2.0


  return d['SC']    
#%%



#%%
#Download data from yahoo finance
start = datetime.datetime(2010,1,1)
end = datetime.datetime(2013,12,31)
ticker = "^GSPC"
d=web.DataReader(ticker,'yahoo',start,end) 

d.info()
#<class 'pandas.core.frame.DataFrame'>
#DatetimeIndex: 1006 entries, 2010-01-04 00:00:00 to 2013-12-31 00:00:00
#Data columns (total 6 columns):
#Open         1006 non-null float64
#High         1006 non-null float64
#Low          1006 non-null float64
#Close        1006 non-null float64
#Volume       1006 non-null int64
#Adj Close    1006 non-null float64


type(d)
#pandas.core.frame.DataFrame


d.head()
d.tail()
#%%



#%%
#calculate KAMA
#---------------
kama = KAMA(d.Close, n=10, pow1=2, pow2=30)

type(kama)
#pandas.core.frame.DataFrame


kama.head(100)
kama.tail(10)
#%%



#%%
df = pd.DataFrame({'price':d.Close, 'KAMA':KAMA(d.Close, n=10, pow1=2, pow2=30) })
df.plot(subplots=True)
#%%

非常感谢。

Answer 1

谢谢你的回复。成功了！

#%%
def KAMA(price, n=10, pow1=2, pow2=30):
    ''' kama indicator '''    
    ''' accepts pandas dataframe of prices '''

    absDiffx = abs(price - price.shift(1) )  

    ER_num = abs( price - price.shift(n) )
    ER_den = pandas.stats.moments.rolling_sum(absDiffx,n)
    ER = ER_num / ER_den

    sc = ( ER*(2.0/(pow1+1)-2.0/(pow2+1.0))+2/(pow2+1.0) ) ** 2.0


    answer = np.zeros(sc.size)
    N = len(answer)
    first_value = True

    for i in range(N):
        if sc[i] != sc[i]:
            answer[i] = np.nan
        else:
            if first_value:
                answer[i] = price[i]
                first_value = False
            else:
                answer[i] = answer[i-1] + sc[i] * (price[i] - answer[i-1])
    return answer
#%%


#%%
#calculate KAMA
#---------------
kama = KAMA(d.Close, n=10, pow1=2, pow2=30)
kama
#%%

Answer 2

您会想要检查这是否真的给出了您所期望的答案，但作为起点，这里有一个非常快速的cython版本。比我更擅长Cython的人可能会发现一些明显的改进。快速比较速度：

In [448]: %timeit df['kama_python'] = filt_python(df['sc'].values, df['price'].values)
10 loops, best of 3: 22.5 ms per loop

In [449]: %timeit df['kama_cython'] = filt_cython(df['sc'].values, df['price'].values)
1000 loops, best of 3: 300 µs per loop

在编写这篇文章时，我从纯python版本开始(见下文)，并不断添加类型注释( pandas docs中的好例子)，直到它变得很快。感谢@Jeff为check np.nan (v != v)指出了正确的方向。

%%cython

cimport numpy as np
import numpy as np
import pandas as pd

cpdef np.ndarray[double] filt_cython(np.ndarray[double] sc, np.ndarray[double] price):
    cdef np.ndarray[double] answer = np.zeros(sc.size)
    cdef int N = len(answer)
    cdef int i
    first_value = True    

    for i in range(N):
        if sc[i] != sc[i]:
            answer[i] = np.nan
        else:
            if first_value:
                answer[i] = price[i]
                first_value = False
            else:
                answer[i] = answer[i-1] + sc[i] * (price[i] - answer[i-1])
    return answer

这是我开始时使用的纯python版本。

def filt_python(sc, price):
    answer = np.zeros(sc.size)
    N = len(answer)
    first_value = True

    for i in range(N):
        if i == 0 or pd.isnull(sc[i]):
            answer[i] = np.nan
        else:
            if first_value:
                answer[i] = price[i]
                first_value = False
            else:
                answer[i] = answer[i-1] + sc[i] * (price[i] - answer[i-1])
    return answer