似たような値動きを探す python

毎日Rといいながら、またpythonです。似たような値動きを探します。

参考にしたのは

https://saturday-in-the-park.netlify.app/TradingTools/09_RecognizeChartShape/

です。感謝です。

　図は、2019-04-25 から 2021-09-09 の 6125岡本工の終値データ から作成したものです。青い太い線は直近30日間のデータです。青太線の最後が 2021-09-09 です。この日付は、凡例の最初の行に示されています。この青太線と同じような値動きをこのデータ内から探すというのがやろうとしていることです。

　凡例の2行目2021-01-20とあります。橙色の細線がその値動きで、2021-01-20に至る30日間の値動きが青太線に近いということです。

　値動きの近さは、相関係数で計ります。2行目r=0.908が相関係数です。この相関係数で近さを計るというのが素晴らしい着想だと思うのです。

　ふつうはそのあとの値動きが気になるので、探し当てた過去データについてはその後の１０日間の値動きを表示しています。

　値動きが近いもの、すなわち相関係数が１に近い方から5日分を描いています。

結構いいところついているように思いますがどうでしょう。

　ソースは次の通りです。

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

def cor(x,y):

  return np.corrcoef(x,y)[0,1]

def similar(df,term=30,after=10,n=5,figsize=(10,10)):

  x=df.reset_index(drop=True)

  v=[(cor(x.iloc[i:i+term].C, x.iloc[-term:].C), i)\

     for i in range(len(x)-term+1-after)]

  v.sort(reverse=True)

  fig,ax=plt.subplots(figsize=figsize)

  ax.plot(x.iloc[-term:].reset_index().C/x.iloc[-1].C,linewidth=3\

          ,label=x.iloc[-1].d)

  for c,i in v[:n]:

    ax.plot(x.iloc[i:i+term+after].reset_index().C/x.iloc[i+term-1].C\

            ,linewidth=1\

            ,label=f"{x.loc[i+term+1].d} r={c:.3f}")

  plt.grid()

  ax.legend()

  return fig,ax

similar(x)

　なお、xはデータフレームで日付を持つd（文字列）と、終値を持つC（数値）が必要です。今回使用したものは次のようなもので、約4年半のデータです。私は株探の時系列データを利用しています。

               d       C

2019-04-25  3035.0

2019-04-26  2965.0

2019-05-07  2890.0

2019-05-08  2987.0

2019-05-09  2868.0

2019-05-10  3175.0

・・・・

2021-09-06  6430.0

2021-09-07  6480.0

2021-09-08  6840.0

2021-09-09  6820.0

　一応プログラムの説明をしてみます。

　相関係数はnumpyのものを使いましたが、そのままでは使いにくいので関数にしました。関数corです。

　similarが主眼の関数です。

　reset_indexが何回か出てきますが、グラフを縦にそろえるため、indexが0からになるように振りなおしています。

　vは1日ずつずらして相関係数をとった結果ですが、相関係数といつからか(i)が必要なので、タプルにして保存しました。この後相関係数でソートするので相関係数を先にしてあります。内包表記です。

　plotのところでは、基準日の値が同じになるように、各データを基準日の値で割っています。これによって、基準日にはすべてのデータが同じ点を通ります。

　以上です。

python

毎日Rといいながら、最近は毎日python生活です。バックテストをしてみたので公開します。
・データは株ドラゴンの時系列データを使いました。マザースの約290銘柄、約１２０日ずつの４本値と出来高です。
・ストーリーは、前日終値より５％下がったら買い、買値より５％上がるか下がったら売るというものです。
・結果はマイナスです。（プラスだったら発表しません）
データは次のようなデータフレームです。列が縦にそろっていなくて見にくくてすみません。右側の５列は自分で付け加えました。

              　　O  　  H    L  　  C     V 　 code    dC     C1      V1 　 V5m  L5H5
d                                                                         
2019-03-20  926  926  908  908   500  3634   NaN    NaN     NaN  NaN   NaN
2019-03-22  907  930  870  905  3400  3634  -3.0  908.0   500.0  NaN   NaN
2019-03-25  950  950  895  896  3500  3634  -9.0  905.0  3400.0  NaN   NaN
2019-03-26  895  915  895  915   900  3634  19.0  896.0  3500.0  NaN   NaN
2019-03-27  960  960  930  956  3800  3634  41.0  915.0   900.0  NaN   NaN

　買の判定は以下の通りです。寄り付きで前日終値の－５%であれば買、寄り付きが買値となります。寄り付きでは前日終値の－５%以下にならなかった場合、日中を見ることになるので、当日安値が前日終値のー５％以下になったら買、買値は前日安値のー５％としました。実際にはこれより安いわけですが。

1. def in_underC(df,k=0.95):
2.   ins = df.C1*k
3.   df["f"]=df.L<ins
4.   df.f=df.f.astype(int)
5.   #in値の設定
6.   df["I"]=-1
7.   df.loc[(df.f)&(df.O>ins),"I"]=ins
8.   df.loc[(df.f)&(df.O<=ins),"I"]=df.O
9.   return df

　売の判定は冗長なのでコードは示さないことにしますが、あらすじは次の通りです。
買値の＋５％を利確値、－５％を損切値ということにします。利確値以上で利確、損切値以下で損切します。
　２日目以降はまず、始値を見て利確、損切になれば始値を売値とします。どちらにも当たらない場合、高値と安値を利確値、損切値と比べます。売りになった場合は利確値、損切値が売値になります。
　買った当日は始値の比較をしないで高値安値との比較のみをします。
　売になった場合は、その値とかかった日数を記録します。買った当日売った場合は０日、売りにならなかった場合マイナス１日とします。

結果のデータフレームです。
            　　f      I 　　 days         g    S　　 state
d                                           
2019-04-25  1  515.0     1 -0.050000  489    ba
2019-04-26  0   -1.0    -1  0.000000    0     -
2019-05-07  1  495.0     0 -0.050000  470    ba
2019-05-08  1  476.0     3  0.050000  499    ba
2019-05-09  1  486.0     2  0.050000  510    ba
2019-05-10  1  478.0     1  0.050000  501    ba
2019-05-13  0   -1.0    -1  0.000000    0     -
2019-05-14  1  483.0     1  0.066253  515     O

先ほどのデータフレームの右側にくっついて出てきますが、その部分だけです。
相変わらず見にくいです。
fは買に入ったかどうかのフラグ
Iは買値。inなのでIを使いました。
daysは売までの日数。入ってないときはー１。０は当日売った場合。
Sは売値
gはS/I-1 いわゆるリターンですが、ゲインという意味でg
stateは売ったのが寄りか場中か。baは場中、最終行はOですが始値で売ったということです。

これを先ほどのデータで検証した結果が次です。

計=36038 f=3732 days=3681 g.mean=-0.0009
100<C<5000 V5m>10000 g.mean=-0.0020
100<C<5000 g.mean=-0.0008
100<C<5000 V5m>10000 5日以内 g.mean=-0.0017
100<C<5000 5日以内 g.mean=-0.0005

計は約２９０銘柄×約１２０日＝36038本のローソク足を検証した。
fは入った日数。約１０分の１入っています。
daysは日数が0以上のものすなわち、売りになった数
g.meanはゲインの平均です。
２行目以降は同じ数字に制限をかけたものです。
100<C<5000は終値が100円より大で5000円より小。100円以下は怪しいし、5000円以上は買えません。V5mは前日までの５日間の出来高平均。１日当たり１００００株くらいの売買がないといけません。
５日以内は５日以内に売れたもののみの集計です。
いずれの場合も売れなかった部分については平均から除外しています。

バックテストでいい数字が出てぬか喜びをしたことは何回もあります。で、心配なので目で見て確認したいです。

ローソク足の上のvは買った日です。

ローソク足中にXがありますが、買った値です。

ローソク足の下にoxがありますが、勝ったか負けたかです。

バックテストだいたい正しいようです。

k-db から hesonogoma で大騒ぎ

　２０１７年末を持ってk-dbが閉鎖ということになりました。お世話になりました。代替サイトですが、hesonogomaさんの株価一覧表を使うこととしました。証券会社で一括データの提供くらいしてくれてもよさそうなものですが、調べた限りでは見つかりませんでした。で、いまだ、対応途上ですが、目鼻が付きましたので、それに伴う様々なトラブルをメモしておきます。

1. データのダウンロード・・・k-dbではurlの指定でcsvがダウンロードできましたが、hesonogomaさんではフラッシュを使用しているようで、そのようなわけにはいかず、手作業でファイルをダウンロードすることとしました。
2. データファイルが４つ必要・・・全銘柄・ETF等・REIT・ファンドと４種類のファイルが必要です。
3. 文字コードの問題１・・・文字コードの問題は相当重症です。ダウンロードしたファイルの文字コードは「UTF-16LE BOM付き」でした。UTF-16をRで読み込むとファイルが途中で切れます。当初原因がわからず、試行錯誤しましたが、StackOverflowに中国人の方が、同じ現象に会って質問をしていました。特定のコード以降読み込まないという趣旨だろうと思います。回答はありませんでした。で、nkf.exeを使っていったんUTF-8に変換することにしました。
4. 文字コードの問題２・・・最初「UTF-8 BOMなし」に変換したのですが、EXCELに読み込むと文字化けする。で、BOM付きに変換。そのせいか、そのファイルをRで読み込むと列数が多いといわれて読み込めない。

ということで現在のところの文字コード部分のコードは次の通り。

system(paste(nkf,"-w8 --overwrite",fnams[1]))
system(paste(nkf,"-w8 --overwrite",fnams[2]))
system(paste(nkf,"-w8 --overwrite",fnams[3]))
system(paste(nkf,"-w8 --overwrite",fnams[4]))
d1=read.csv(fnams[1],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d2=read.csv(fnams[2],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d3=read.csv(fnams[3],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d4=read.csv(fnams[4],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)

なお、「UTF-8 BOM付き」をEXCELで読み込むと、tabで区切ってくれないのでもうひと手間かかります。

また、fileEncodingに「UTF-8-BOM」はありますが「UTF-16-BOM」はありませんでした。「UTF-16」はありますが上記の問題に会います。

さらに、４つのファイルの項目名が違う。全銘柄ファイルでは業種、それ以外のファイルでは種別となっている。

その部分を乗り越え、一つのファイルにするのが次のコード

names(d2)[which(names(d2)=="種別")]="業種" #names(d2)=gsub("種別","業種",names(d2))のほうがおしゃれかな？
names(d3)[which(names(d3)=="種別")]="業種"
names(d4)[which(names(d4)=="種別")]="業種"
orinames=c("SC","名称","市場","業種","始値","高値","安値","株価","出来高","売買代金.千円.")
d=rbind(d1[,orinames],d2[,orinames],d3[,orinames],d4[,orinames])

とりあえず、一つのデータフレームにできました。次の問題は

1. 銘柄名が長い。・・・新しい銘柄名に変えようか、迷っている。
2. 市場名が違う。・・・これは古いものに統一しようと思っている。
3. 銘柄コードが違う。・・・k-dbでは「1301-T」のような形式でしたが新しいほうは「1301」のような形式です。新しい形式にしたときに問題になりそうなことは、
• データ参照の時にyt[["1301-T"]]という参照方法とyt[[1]]という参照方法がありましたが、「1301」形式では文字列として扱われるのか整数として扱われるのか、悩ましい。Rが勝手に変換してくれるのもよいのやら。
• k-dbでは同一コード同一銘柄で別市場のものがありましたが、新しいほうでは同一コードが見当たりません。東証と別市場がある場合は東証のみになっているのでしょうか？

ということで、大晦日から元旦にかけての作業でした。

R ディレクトリのみ表示

ディレクトリのみ表示する。

t=dir()
t[file.info(t)$isdir]

１文にするなら
dir()[file.info(dir())$isdir]

こんな命令があった。
list.dirs(".",recursive = F)

参考https://stackoverflow.com/questions/4749783/how-to-obtain-a-list-of-directories-within-a-directory-like-list-files-but-i

R 関数を文字列で呼びだす

do.call("mean",list(1:10,na.rm=T))
第２引数は必ずリストで。これが関数への引数リストとなる。

eval((parse(text= "mean(1:10,na.rm=T)")))
関数に限らず文字列を実行したいとき。

R 練習問題

Rの練習問題を作ってみました。今回の対象は次のようなデータフレームです。試験の点数処理を題材としてみました。データを作るスクリプトです。乱数でとも思いましたが、同じ結果の方が見やすいでしょうから。

データはこんな風です。３組それぞれ５人の国・数・英の点数です。

問題と出力です。
#問１　dに個人ごとの計を付加

#問２　計順にソートして表示

#問３　dに計で学年順位を付加

#問４　組ごと計順に表示

#問５　教科間の相関係数

#問６　科目別箱髭図

#問７　国の組別箱髭図

#問８　教科集計

#問９　計の度数分布

#問１０　科目別平均

#問１１　計の組別平均

#問１２　科目別組別平均

#問１３　国の組別平均・最高・最低

#問１４　国50未満を抽出

#問１５　国50未満を抽出、組番国のみ表示

#問１６　数40以上英50以上を抽出

#問１７　数30以上を計順に表示

#問１８　科目別、成績順、組番

#問１９　計の２０点刻み度数分布図

#問２０　科目別５点刻み度数分布表


解答


終わりに
この練習問題では科目は横に展開し、組は縦に展開しているので、科目ごとと組ごとの扱いが違って来る。
apply系の命令をうまく使うことは難しい。この練習問題では、結構面白い使い方ができたと自負している。

簡単にポイントにふれておく。

問１　①新たな列の付けくわえ方②d[3:5]はデータフレームdの３列目から５列目に限定している。データフレームはベクトルのリストであるので、d[3:5]で３列目から５列目が取り出される。これがmatrixならば、d[,3:5]とカンマとつけなければならない。もちろんデータフレームでもこのように指定してもよいが。③行集計にrowSumsを用いた。apply(d[3:5],1,sum)としてもよい。

問２　①ソートにはorderを用いるのが良い。単一のベクトルをソートするにはsortを用いればよいが、キーでソートするようなときはorderがよい。orderは昇順にソートするので、高得点の方が先に来るように、マイナスした。desc=Tを用いてもよいが、おしゃれじゃない。②行番号で指定するので、d[order(-d$計),]と、コンマが必要となる。問１で述べたように、ここでコンマがないと、列の指定となってしまう。

問３　①順位はrank。②引数のt=は ties.method = c("average", "first", "random", "max", "min")である。引数の変数名は誤解のない範囲で省略できる。

問４　orderでソート項目は並べればよい。

問５　特にない

問６　これだけで科目ごとの箱髭図を書いてくれるなんて、なんて便利なんだろう。

問７　組ごとの箱髭図はこう指定する。問６との違いは、先に述べたように、科目は横に展開し、組は縦に展開していること。

問８　特にない

問９　特にない

問１０　d[3:5]がデータフレームなので、sapplyはリストの要素の各ベクトルすなわち各列に対して動作する。

問１１　問１０との違いは、前述のとおり

問１２　applyの中でtapplyを使ってみた。うまくいった。sapplyでもいいのだね。

問１３　いろいろと試行錯誤し、もっともよさげなものをあげた。

問１４～１７　抽出

問１７　抽出と整列を同時にできないかというチャレンジ。一時変数を使わないでやりたかったが、最終的には用いてしまった。敗北感あり。使わない版もいくつか作ったが、同じ計算を２度やらせていて、許せない。いつの日か思いつくかもしれない。

問１８　科目別のベストテンなどの表はよく作るだろう。それである。文字列ベクトルを作っておいて、その並び順を添え字で指定するというこのパターンは初めての経験かもしれない。

問１９　breaksについては問２０参照

問２０　科目別度数分布表。苦労したのは、区間の指定。たとえば、breaks=c(10,20,30,40)とすると、この４つの数値によって、３つの区間に分けられる。10≦x≦20,20＜x≦30,30＜x≦40。最初の10の不等号に＝が入っている

サンプルデータでは、回帰係数や回帰分析や、主成分分析、クラスタ分析も面白くなかろうと思ってやめた。

移動平均について

移動平均について考察してみる。なお、単純な移動平均についてである。
本日の終値を C0、前日の終値をC1、前々日の終値をC2、…と書く。
騰落率は C0/C1-1である。
５日移動平均MA5=(C0+C1+C2+C3+C4)/5である。

移動平均からの乖離率について

今、１日の騰落率が５日間ともrであったとしよう。すると、
C0 = C4×(1+r)^4 ≒ C4×(1+4r)
ｒが十分0に近いとき、(1+r)^n ≒ 1+nr であることを用いている。 今後これをどしどし使う。

同様に、C1≒C4×(1+3r)、C2≒C4×(1+2r)、C3≒C4×(1+r)、であるから、
MA5 ≒ C4×{(1+4r)+(1+3r)+(1+2r)+(1+r)+1}/5 ≒ C4×(1+2r) となる。
C0のMA5からの乖離率は、
C0/MA5-1 ≒ (1+4r)/(1+2r)-1
　　　　　　　≒ (1+4r)(1-2r)-1
　　　　　　　≒ 2r
すなわち、騰落率の２倍の値となる。逆に言えば、乖離率の２分の１が１日当たりの騰落率ということになる。
同様に、（２５日移動平均からの乖離率）÷１２ ＝ ２５日間の平均騰落率
　　　　　（７５日移動平均からの乖離率）÷３７ ＝ ７５日間の平均騰落率
と考えることができる。

移動平均曲線の傾き

移動平均が右肩上がりならば、上昇トレンドにあると考えられる。しかし、チャートを比べて上昇の度合いを考えたいと思っても、チャートの縦の比率が一定ではないので、単純には比べられない。
移動平均曲線の傾きはどう求めたらよいだろうか。前述のMA5で調べてみる。

最終のMA5の傾きは、前日MA5と当日MA5の差である。
(MA5の傾き)＝(当日MA5)-(前日MA5)
　　　　　　　　＝ (C0+C1+C2+C3+C4)/5-(C1+C2+C3+C4+C5)/5
　　　　　　　　＝ (C0-C5)/5

 なんと、５日前との終値の差の５分の１である。途中の値は関係ないのだ。
C0-C5はモメンタムというらしい。
同様に、２５日移動平均の最終日の傾きは、２５日前との終値の差の２５分の１となる。
ただの終値の差に、こんな意味があったとは！！

さて、実際には、株価が高いところでは当然その傾きも大きくなるのだから、比較のためには、
C5/C0-1の値を使うことになるであろう。この値はROC(rate of change)というらしい。

quantmodでは momentum(x, n = 1, na.pad = TRUE)
ROC(x, n = 1, type = c("continuous", "discrete", na.pad = TRUE)
となっている。上の計算と合わせるなら"discrete"である。