前回、前々回とシストレ入門として、STC（Stocastic Oscilator）について、記載してきた。
今回は、さらに精度が高い(売買シグナルがより速く出る)というMACD（Moving Average Convergence Divergence）をPythonで再現し、さらにDecomposeすることにより分かり易い判断が出来るようしたいと思う。以下のようなグラフが得られます。
※横軸を日付にしたグラフを示します
以下はDecomposeしてTrendのみを利用して元の変化と一緒に描画したものです。非常にシンプルになったので、売買のシグナルが分かり易いです

以下は、普通のMACDグラフです。

【参考】
①ＭＡＣＤとは？売買シグナルが早く出現する、移動平均線の進化系!?
②Pythonで書いた指数移動平均(EMA)のコードの比較
③移動平均@wikipedia
###やったこと
・StocasticOscilatorの進化系にしてみた
・MACDの計算に出てくるEMAについて
・MACDのコード解説
###・StocasticOscilatorの進化系にしてみた
前回は、Decompose導入して見やすいとはいえ、以下のグラフから読み取って判断していた。

ここでは、20%以下は売られすぎ、80%以上は買われすぎという判断をしつつ、上昇場面では青線が上に来た変化点が買いサイン、逆にオレンジ線が上に来た時が売りサインと判断する。
これをシステマティックに判断するために、この20%、80%の要素を入れるために、以下の関数で(%K-%D)/%Kを計算して判断できるようにした。
※因子１は、0devideを回避するために入れている

def calc_dff(dK,dD):
    dff= dK.copy()
    for j in range(len(dK)):
        #print(j)
        if dK[j]-dD[j]>0:
            dff[j]=100*(dK[j]-dD[j])/(abs(dK[j])+1)
        else:
            dff[j]=100*(dK[j]-dD[j])/(100-abs(dK[j])+1)
    return dff

この関数で計算すると以下の一番下のグラフが得られる。
このグラフを見ると、シグナルはより抑制され売り時買い時が明確になっている。これが通常の単純移動平均の限界だと思われる。

したがって、上記のグラフは以下を示している。

Trend＝Observed－Cycle
5日SMA（単純移動平均）、25日SMA
売買シグナル＝100*(5日SMA-25日SMA)/5日SMA

###・MACDの計算に出てくるEMAについて
EMA( Exponential Moving Average)について、書いておこう。
定義は参考③にあるように以下のとおり
t≧3 の場合の EMA の計算式は次のとおりである。

S_{{t}}=\alpha \times Y_{{t-1}}+(1-\alpha )\times S_{{t-1}}

ここで、$S_{{t}}$はある時点 t での EMA であり、$Y_{{t-1}}$は時系列上のある時点 t の観測値である。また、近似的に$\alpha = {2 \over {N+1}}$である。
このEMAを利用して、MACD等は以下のように定義される。

MACD＝12日EMA－26日EMA
Signal＝MACDの9日EMA
ヒストグラム＝Signal－MACD

したがって、上記のStocasticOscilatorと比較すると、単純移動平均を指数移動平均に置き換えて長期と短期移動平均の日数を少し変更したものと考えることが出来る。
ヒストグラムが上記の一番下のグラフに相当している。
###・MACDのコード解説
以下のコードが今回のシストレの最終結果だと思います。
参考②のコードを利用した分かり易いものを説明します。
※横軸を日付にした最初のグラフを描画するコードはおまけに掲載します
利用するLibは以下のとおりです。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import datetime as dt
from pandas_datareader import data
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL

以下はデータ取得の関数です。

def get_stock(stock,start,end):
    df = data.DataReader(stock, 'stooq',start)["Close"]
    df = df.iloc[::-1]
    return df[start:end]

パラメーター設定します。

stock0 = '9437'
stock = stock0 + '.JP'
bunseki = "series"
start = dt.date(2020,1,1)
end = dt.date(2020,6,5)
df = pd.DataFrame(get_stock(stock, start, end))
series=df['Close']
cycle, trend = sm.tsa.filters.hpfilter(series, 144)
df['Close']=trend

df['Close'].values.tolist()でlistにデータを格納します。

gwalk = df['Close'].values.tolist()
print(gwalk[0:3])

[2978.67, 3016.01, 2972.77]のようなデータが出力されます。
以下でn日EMAを計算します。参考②ではaを0.15に固定していましたが、ここではaを計算することが大切です。

def EMA1(x, n):
    a= 2/(n+1)
    return pd.Series(x).ewm(alpha=a).mean()

上記の関数により、12日EMAなどは簡単に計算できます。
また、MACDやSignalも以下のように簡単に計算できます。

y12 = EMA1(gwalk, 12)
y26 = EMA1(gwalk, 26)
MACD = y12 -y26
signal = EMA1(MACD, 9)
hist_=MACD-signal

ここでは横軸を数値列で表現します。明示的にplotの横軸として記載してもいいですが、ここでは不要です。
後はグラフを描くだけです。
とりあえず、ax2.bar(ind,hist_)が面倒なところです。

ind = np.arange(len(signal))
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 12, 4*2),dpi=200)
ax1.plot(gwalk,label="series")
ax1.plot(y12,label="y12")
ax1.plot(y26,label="y26")
ax2.plot(MACD,label="MACD")
ax2.plot(signal,label="signal")
ax2.bar(ind,hist_)
ax1.legend()
ax2.legend()
ax1.grid()
ax2.grid()
plt.show()
plt.savefig("ema_.png")
plt.close()

この結果は以下のように出力されます。

decomposeした結果は以下のようにシンプルです。
以下は9437NTTドコモの結果ですが、もうすぐプラテンしそうです。

つまり、ノイズが乗ろうが無かろうが、要はヒストグラムがマイナスからプラスに変化する部分で買いだし、逆にプラスからマイナスに変化する部分が売りだということで、ほぼ実際の動きにあっているように見えます。
ロングスパンの売買には有効な指標だということが分かります。
もう一つの横軸を日付で描くものをおまけで説明します。
###まとめ
・StocasticOscilatorを分かり易い指標で表現してみた
・MACDを利用した指標はロングスパン売買に有効そうである

・やはり、ここまで分かり易いのでいわゆるシストレに応用したい
###おまけ
あまりコード的には変わりませんが、横軸を日付にするのはとても苦労したので、掲載しておきます。じっくり読めば楽しめると思います。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
from pandas_datareader import data
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL

def get_stock(stock,start,end):
    df = data.DataReader(stock, 'stooq',start)["Close"]
    df = df.iloc[::-1]
    return df[start:end]

from numba import jit
@jit(nopython=True)
def EMA3(x, n):
    alpha = 2/(n+1)
    y = np.empty_like(x)
    y[0] = x[0]
    for i in range(1,len(x)):
        y[i] = alpha*x[i] + (1-alpha)*y[i-1]
    return y

def EMA1(x, n):
    a= 2/(n+1)
    return pd.Series(x).ewm(alpha=a).mean()

def calc_dff(dK,dD):
    dff= dK.copy()
    for j in range(len(dK)):
        dff[j]=100*(dK[j]-dD[j])/(abs(dK[j])+abs(dD[j])+1)
    return dff

stock0 = '9437'
stock = stock0 + '.JP'
start = dt.date(2020,1,1)
end = dt.date(2020,6,5)
df = pd.DataFrame(get_stock(stock, start, end))
date_df=df['Close'].index.tolist() #ここがポイント
print(date_df[0:30])
series = df['Close'].values.tolist()

bunseki = "trend" #series" #cycle" #trend
cycle, trend = sm.tsa.filters.hpfilter(series, 144)
series2 = trend

y12 = EMA3(series2, 12)
y26 = EMA3(series2, 26)
MACD = y12 -y26
signal = EMA3(MACD, 9)
hist_=MACD-signal

ind3=date_df[:]
print(len(series),len(ind3))

fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 8, 4*2),dpi=200)
ax1.plot(ind3,series,label="series")
ax1.plot(ind3,series2,label="series2")
ax1.plot(ind3,y12,label="y12")
ax1.plot(ind3,y26,label="y26")
ax2.plot(ind3,MACD,label="MACD")
ax2.plot(ind3,signal,label="signal")
ax2.bar(ind3,hist_)
ax1.legend()
ax2.legend()
ax1.grid()
ax2.grid()
plt.savefig("./stock/{}/ema_decompose_%5K%25D_{}_{}now{}.png".format(stock0,stock,bunseki,start))
plt.pause(1)
plt.close()

df['Close']=series  #series" #cycle" #trend
df['series2']=series2
df['y12'] = EMA1(df['Close'], 12)
df['y26'] =  EMA1(df['Close'], 26)
df['MACD'] = df['y12'] -df['y26']
df['signal'] = EMA1(df['MACD'], 9)
df['hist_']=df['MACD']-df['signal']
date_df=df['Close'].index.tolist()
print(df[0:30])

fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 8, 4*2),dpi=200)
ax1.plot(df['Close'],label="series")
ax1.plot(df['series2'],label="series2")
ax1.plot(df['y12'],label="y12")
ax1.plot(df['y26'],label="y26")
ax2.plot(df['MACD'],label="MACD")
ax2.plot(df['signal'],label="signal")
ax2.bar(date_df,df['hist_'])
ax1.legend()
ax2.legend()
ax1.grid()
ax2.grid()
plt.savefig("./stock/{}/ema_df_decompose_%5K%25D_{}_{}now{}.png".format(stock0,stock,bunseki,start))
plt.pause(1)
plt.close()