UNIX timestamp
まとめ
本文
UNIX timestampをちょっとだけまとめ。
PythonでUNIX時間(エポック秒)と日時datetimeを相互変換 | note.nkmk.me
datetime.datetime.fromtimestamp
は実行環境によってtimezoneが変わるみたいなので、timezoneは指定しておいた方が無難そう。どのtimezoneに指定するかは触っているデータによって変わる気が。日本人ユーザーの消費行動ならJSTに合わせた方がいいし、外国ならその地の時間に合わせた方がよさそう。今回私が回しているkaggleの環境設定はUTCでした。
print(datetime.datetime.fromtimestamp(0)) # machineの標準設定 print(datetime.datetime.fromtimestamp(0, datetime.timezone.utc)) # timezone 指定 print(datetime.datetime.fromtimestamp(0, datetime.timezone(datetime.timedelta(hours= 9))))
変換自体は余裕だと思ったらミリ秒まで含む含まないで10桁or13桁の2種類あるみたいです。これはフェイントですね。
timestamp(13桁)カラムをDateTime型に変換したい
- 13桁のまま突っ込む場合(エラー)
%%time a = train.iloc[0, 2] print(a) print(datetime.datetime.fromtimestamp(a, datetime.timezone.utc))
1659304800025
ValueError Traceback (most recent call last)
ValueError: year 54551 is out of range
- 10桁として突っ込む場合(OK)
%%time unix_adjust = 1000 # ミリ秒まではいらないので a = train.iloc[0, 2] print(a) print(datetime.datetime.fromtimestamp(a / unix_adjust, datetime.timezone.utc))
1659304800025 2022-07-31 22:00:00.025000+00:00 CPU times: user 7.86 ms, sys: 989 µs, total: 8.85 ms Wall time: 8.02 ms
せっかくなのでUNIX timestampをdatetimeに変換するところまでやってみます。データ数が22mio超えと大きいのでapply
のnp.vectorize
を使用したときの速度の違いもみてみみます。5秒くらいnp.vectorize
を使用することで短縮できました、データ数が増えればもっと短くなるんですかね?
apply
を使用したとき
%%time print(len(train)) unix_adjust = 1000 # ミリ秒まではいらないので # apply train["ts_datetime"] = train["ts"].apply(lambda x: datetime.datetime.fromtimestamp(x / unix_adjust, datetime.timezone.utc))
22547356 CPU times: user 30 s, sys: 1.73 s, total: 31.7 s Wall time: 31.6 s
np.vectorize
を使用したとき
%%time print(len(train)) unix_adjust = 1000 # ミリ秒まではいらないので # np.vectorize def func_0(ts): return datetime.datetime.fromtimestamp(ts / unix_adjust, datetime.timezone.utc) train["ts_datetime"] = np.vectorize(func_0)(train["ts"])
22547356 CPU times: user 23.4 s, sys: 3.38 s, total: 26.8 s Wall time: 26.8 s
jsonlの読み込み
まとめ
- jsonlファイルで中身がバカでかいやつは一気に読み込むとメモリエラーとなるのでchunkを使用するとよい。
- jsonの中身は辞書型みたいな感じでデータ格納されてるっぽい。
- バカでかいファイルの読み込みは一度parquet形式で保存して次回からそれを読み込むと速い。
本文
jsonもあんまりよくわかってませんが、jsonlが出てきたのでメモ。
pandasでJSON文字列・ファイルを読み込み(read_json) | note.nkmk.me
まずは上記サイト上の「JSON Lines(.jsonl)を読み込み」を参考にファイルを読み込んでみましょう。
%%time test = pd.read_json(DATA / "test.jsonl", orient= "record", lines= True) print(type(test)) display(test.head())
dataframe形式で返ってきました。ただこれは読み込む先のjsonlファイル上の問題ですが、events列が以下の通りlistとなっていてさらにその中身がdictとなっていて扱いずらそうな。
print(test.iloc[0,1]) print(type(test.iloc[0,1])) print(type(test.iloc[0,1][0]))
そこで、他の人がやっていたchunkを使用して読み込むパターンをやってみましょう。ドキュメントを見る感じchuncksize
はjsonlを読み込むときのみに使用できるみたい。ただ何を表しているのかいまいちわからなかったので、手を動かしながらやってみます。
pandas.read_json — pandas 1.5.1 documentation
%%time cnt = 0 CHUNKSIZE = 10000 test = pd.read_json(DATA / "test.jsonl", lines= True, chunksize= CHUNKSIZE) print(type(test)) for chunk in test: cnt += 1 if cnt > 1: continue display(chunk) print(type(chunk)) print(cnt)
うーん、まず返してくるのはドキュメントの通りJsonReaderという形式です。for文とかで一行ずつ取り出してみるとchuncksize
を使用しなかったときのようなdataframe形式で取得できました。1chunkあたりのデータ数がちょうどCHUNKSIZE
と一致しているので、chucksize
は1回?あたりに読み込む行数みたいなイメージでしょうか。また、今回の場合はfor文168回回しているので、全データ数は10000*168=1,680,000行くらいでしょうか。大きいですね。
結局なぜこのchuncksize
を使用してJsonReader形式のデータを返しているのかとう話ですが、おそらく返ってくる値がdictとかでめんどくさいのを読み込みの段階でうまく扱えるからですかね?読み込む早さはchuncksize
使用している方が読み込んだ後にfor文回しているので遅かったです(ちなむと読み込み自体は以下の通りJsonReader形式で読み込んだ方が速かったです。)そもそもjsonオブジェクトって項目名と値を:で結んだりして表現するみたいなので、辞書型みたいで複雑になるのはあるあるなのでしょうか。また、それらをそのままの形式で一気に抽出してからいじるより、サクッとjsonreader形式で抽出してからfor文で辞書型の中身をいじった方が早かったりするんですかね?
%%time CHUNKSIZE = 10000 test = pd.read_json(DATA / "test.jsonl", lines= True, chunksize= CHUNKSIZE)
jsonreader + for文でaid
項目のみ抽出したとき(30s)
%%time CHUNKSIZE = 10000 test = pd.read_json(DATA / "test.jsonl", lines= True, chunksize= CHUNKSIZE) print(type(test)) for chunk in test: for session_id, events in chunk.values: last_aid_list.append([session_id, events[-1]["aid"]])
lines=Trueで抽出してからapplyで加工したとき(23s)
%%time test = pd.read_json(DATA / "test.jsonl", orient= "record", lines= True) test["aid"] = test["events"].apply(lambda x: x[-1]["aid"])
時間測ってみたもののむしろapplyで加工した方が速かった。。うーん、jsonreaderで読み込むことにどういうメリットがあるんだろう。。
上記までもそれなりのデータ数がでしたが、さらに大きいデータになると一気にloadしようとする方はメモリエラー出た一方でchunck
は読み込めました。メモリ削減の効果はありそうです。ただ、これだけのデータ数になると読み込みにとにかく時間がかかる。。8分読み込んでもまだ。。
%%time def read_data(path, is_train): ''' data読み込み ''' # # lines 読み込み # # データ数大きいとメモリ不足 # df = pd.read_json(DATA / path, orient= "record", lines= True) # df["aid"] = df["events"].apply(lambda x: x[-1]["aid"]) # df["ts"] = df["events"].apply(lambda x: x[-1]["ts"]) # df["type"] = df["events"].apply(lambda x: x[-1]["type"]) # jsonreader 読み込み last_aid_list = [] CHUNKSIZE = 10000 json_reader = pd.read_json(DATA / path, lines= True, chunksize= CHUNKSIZE) for chunk in json_reader: for session_id, events in chunk.values: last_aid_list.append([session_id, events[-1]["aid"]]) # del json_reader,chunk del chunk gc.collect() df = pd.DataFrame(last_aid_list) del last_aid_list gc.collect() return df train = read_data("train.jsonl", is_train= True)
ということで、parquet formatで早く読み込むというnotebookがあったのでこちらを使わせて頂く。parquetに一度保存して次回からはparquetを読み込むという方法。確かにこれなら読み込みにかかる時間がだいぶ減った。にしても5 / 124ファイルのデータ数が22,547,356行とは。。ひいい。
%%time def read_data(path, is_train, is_save= True): ''' data読み込み ''' # =================== # read existing parquet files files = sorted(glob('../input/otto-chunk-data-inparquet-format/train_parquet/*'))[0:5] # files = sorted(glob('../input/otto-chunk-data-inparquet-format/train_parquet/*')) dfs = [] for path in tqdm(files): dfs.append(pd.read_parquet(path)) dfs = pd.concat(dfs).reset_index(drop=True) return dfs train = read_data("train.jsonl", is_train= True)
おしまい。
連続型分布と標本分布
自分用メモです。間違えていることあります。
〇指数分布とガンマ分布
↑体重がガンマ分布に従わない?従う従わないの判定方法は?
身長は正規分布、体重はガンマ分布…なぜ? - ほうかいのじゅもん
↑「正の値をとる確率変数はガンマ分布に従う」?shape大→正規分布へ。
↑「体重=縦横高さの積で近似」「ガンマ分布の3乗根は正規分布」「身長が正規分布」→「体重はガンマ分布」?
〇多変量正規分布
↑導出
2変量正規分布の条件付き確率の期待値と分散(統計検定準1級対策3)
条件付き分布
重回帰分析参考サイト
自分用のメモです。言葉足らずですし、間違えていることもあります。
重回帰モデルの理論と実装 -なぜ正則化が必要か- | Deploy on Friday
正則化ってなにをしてるかについて。正規方程式でが正則でない(逆行列もたない)ときになんらかで調整したい。例えばリッジは誤差二乗和にを足すことで、で微分したりして進めたときの正規方程式の左辺がみたいな感じになる。
幾何的にみた重回帰の解釈。にを作用させたら勝手にに射影してくれて、射影した先 - とは直交してくれるものだと思ってた。あまりまだ正しく理解できていない気がする。
正規方程式(normal equation)の導出と直感的理解 - あつまれ統計の森
正規方程式の導出がよくわからないときにみた。ベクトルによる微分の説明もあり。
27-6. 回帰の有意性の検定 | 統計学の時間 | 統計WEB
回帰式自体に意味があるかの検定。言っていることはなんとなくわかるのだが、若干まだもやっと。なんとなくは。。
自己相関の検定
時系列データが自己相関を持つか=標本自己相関≠0の検定。与えられた時系列データがiid系列と仮定すれば標本自己相関の漸近分布は平均0分散1/Tの正規分布に従う。iid系列という仮定が結構強い気もするけど、とりあえずやってみる。
参考サイト
statsmodels.graphics.tsaplots.plot_acf — statsmodels
Pythonで自己相関グラフ(コレログラム)を描く - Qiita
やってみた
pd.plottingを使うなどの方法もあったけど、パラメータ設定がやりやすそうだったので今回はstatsmodelsで行うことに。それっぽいものはできたけど、95%信頼区間はBartlett’s formulaなるものを用いているらしい。ARMA過程やらの考えを用いているみたいだけど、ちょっとドツボにはまりそうだったので、おいおいということで一旦はそれっぽいものとしておいておく。。でも、ネットで落ちているようなコレログラムをみると信頼区間はラグによって結構動いているけど、今回はどれもほぼほぼ一定に見えるのでもしかしたらこれは特徴的?
from statsmodels.graphics.tsaplots import plot_acf name= "USDJPY Curncy" col = "log_chg" tmp = df[df.id == name]\ .sort_values("date")\ .loc[:, col] plot_acf( x= tmp.dropna().values, lags= 20, # ラグいくつまで表示するか alpha= 0.05, # 信頼区間 zero= False, # lag= 0を表示するか title= name.split(" ")[0] + " Autocorrelation", # title # missing= "drop", # トキュメントにはあるのに使えなかった );
例の通り、各資産についてやってみた。自己相関がそんなに強いものはあまりない印象。TPXやUSDJPYは逆相関強いという噂はなんとなく聞いていたけど、過去20年に限ってはTPXがラグ1-3でそれなりにあるかな?という感じ。SP500(ES)のラグ1が最も大きくしかも逆相関なのは意外。大きく下落する前に極まって最後の上昇を見せることがあるけど、その影響が強く出ている?米債(TY)のラグ20も比較的大きく逆相関なのは、20って1カ月くらいだからちょっと気になる。
name = "log_chg" fig, ax = plt.subplots(3,2, figsize= (15,12)) for i, id_ in enumerate(param['ticker']): tmp = df[df.id == id_]\ .sort_values("date")\ .dropna()\ .loc[:, name]\ .values plot_acf( x= tmp, lags= 20, # ラグいくつまで表示するか alpha= 0.05, # 信頼区間 zero= False, # lag= 0を表示するか title= id_.split(" ")[0] + " Autocorrelation", # title ax = ax[i // 2][i % 2], ); ax[i // 2][i % 2].set_ylim(-.12, .12) plt.show() # fig.savefig("figs/autcorr.png")
SPとTPのラグ0とラグ1を散布図でプロットしてみると、5%以上動いた際の挙動で自己相関が大きく決定している気が。なので、ラグ0,ラグ11ともに5%以下のときに絞ってもコレログラムを作成してみた。5%以上動かないときが99.5%なのでほぼこの範囲内に含まれている。すると、それほど特徴のない結果に。ほとんど有意に自己相関を持つ資産やラグがない。。3%や2%も試したけど似たような結果でした。
id_ = "ES1 Index" col = "log_chg" # 自己相関算出 shift= 1 series = df[df.id == id_]\ .sort_values("date")\ .loc[:, col]\ .dropna() corr = np.corrcoef( series.values[shift:], np.roll(series.values, shift= shift)[shift:], ) print(corr) # ラグnでプロット fig, ax = plt.subplots() pd.plotting.lag_plot(series, lag=shift, ax= ax) ax.set_xlim([-.15, .15]) ax.set_ylim([-.15, .15]) ax.grid(axis= "x") ax.grid(axis= "y") ax.set_title("SPX") fig.show()
TPX版
df_small = df[ (df.log_chg.abs() <= .05) & (df.log_chg_lag_1.abs() <= .05) ] print(len(df.dropna(subset= ["log_chg", "log_chg_lag_1"]))) print(len(df_small)) print(len(df_small) / len(df.dropna(subset= ["log_chg", "log_chg_lag_1"]))) name = "log_chg" fig, ax = plt.subplots(3,2, figsize= (15,12)) for i, id_ in enumerate(param['ticker']): # tmp = df[df.id == id_]\ # .sort_values("date")\ # .dropna()\ # .loc[:, name]\ # .values tmp = df_small[df_small.id == id_]\ .sort_values("date")\ .dropna()\ .loc[:, name]\ .values plot_acf( x= tmp, lags= 20, # ラグいくつまで表示するか alpha= 0.05, # 信頼区間 zero= False, # lag= 0を表示するか title= id_.split(" ")[0] + " Autocorrelation", # title ax = ax[i // 2][i % 2], ); ax[i // 2][i % 2].set_ylim(-.12, .12) plt.show() # fig.savefig("figs/autcorr.png")
まとめ
株価等のデータは裾が厚いという特徴あり+相関は外れ値の影響受けやすい=全データを対象に自己相関について判断すると米株は逆相関が強いみたいなミスリードが生じるので注意が必要ですな。
特に株価が3%ほど大きく動いた際に米株が逆相関っぽくなっているのは、おそらく米株が下落後に反転する力が強いとかが原因な気が。半面、コロナやリーマンのときらへんに前日大きく上昇&翌日大きく下落のような場面もみるので、まあこんなところが逆相関っぽくみえるところに影響したのではないかと思います。
Q-Qプロット
株やら為替のリターンは、正規分布より裾が厚い分布になっているとなんとなく教えられてきたけど、本当にそうなの?と思って調べてみた。与えられた分布の正規性の確認方法は1.ヒストグラム2.Q-Qプロット3.検定(shapiro-wilk or kolmogorov-smirnov)のどれかみたいだけど、今回は2.Q-Qプロットをやってみる。
参考サイト
Q-Qプロットを理解するためPythonで描いてみる - Qiita
statsmodels.graphics.gofplots.qqplot — statsmodels
【Python】QQプロットとは?対数変換とは?|初心者向けに解説 | 月見ブログ
やってみた
statsmodels.api
内のqqplot
使うのが一番簡単そう。なので、使ったうえで、パラメータの挙動諸々を確認してみる。
使用データ
ドル円の対数リターンを使います。
実際のコード
Q-Qプロットなのでやっているのは「与えられた分布」と「設定した分布(正規分布など)」の分布が一致していそうかの確認。そのために、各分布上の点の分位数を算出し、一致するもの同士をプロットしていく。
とりあえず参考サイト上にあったコードをそのままコピペしてそれっぽく変数を変えてみます。このままでは、「ハア。。??」という感じではあるものの、上下に裾が大きくなっていそうなことは確認できる。
import statsmodels.api as sm import scipy.stats as stats id_ = "USDJPY Curncy" name = "log_chg" fig, ax = plt.subplots(1,1) sm.qqplot( df.loc[df.id == id_, name].dropna(), dist= stats.t, distargs=(4,), # 自由度 fit=True, line="45", ax = ax ) plt.show()
stats.modelのドキュメントを見ながら、パラメータなどについて確認してみます。主だったパラメータは以下の通りだと思います。
不変分散を基に母集団について検定するので、設定する分布はt分布ですね。ただ、サンプル数が5217と十分大きいので、ほぼ正規分布に一致しそうではある。
import statsmodels.api as sm import scipy.stats as stats id_ = "USDJPY Curncy" name = "log_chg" fig, ax = plt.subplots(1,1) sm.qqplot( data= df.loc[df.id == id_, name].dropna(), dist= stats.t, # 比較する分布 distargs=(len(df.loc[df.id == id_, name].dropna() - 1),), # distで設定した分布へのパラメータ(今回は自由度) fit=False, # dataをスケール line="s", ax = ax ) plt.show()
最後に作成したプロットをstatsmodel
を使わずに再現してみることで何しているのか確認してみる。いい感じに再現できているでしょう。サンプル数が5000くらい。
# Q-Qプロット手作成 # Q-Qプロット手作成 tmp = df.loc[df.id == id_, name] n=len(tmp) m=n+1 ii=np.arange(n)+1 # i=1~n qo=np.sort(tmp) qp=stats.norm.ppf( ii/m ) qpp=stats.norm.rvs( size= len(tmp), loc= tmp.mean(), scale= tmp.std(), ) # plot plt.plot(qp, np.sort(qpp)) plt.plot(qp,qo,'r*') plt.ylabel('Observed Quartile') plt.xlabel('Predicted Quartile') # # plt.savefig('qqplot.png',dpi=300)
こうみると、上も(正規分布に属しているなら)もっと小さな値で上位数パーセントの数値が出てくるところを、より大きな値で上位数パーセントが出てしまっている=右側の裾があつい、と解釈。下側も同様に解釈して左右の裾が正規分布と比較してあつくなっていることをわかりやすく表せました。
ドル円以外でもQ-Qプロットを作成してみる
ドル円以外はどうなのでしょう??TPXは下落の方が裾が大きくなっているみたいだけど、SPXはそんなことなさそう。(下落時の方が値幅は大きくなりそうなので、SPXの結果は意外)。為替はドル円は左右でそんなに変わらないけどEURUSDは上昇の方が大きいのが意外。リスクオフではドル買いでは。債券は米債欧州債ともに売られる方が値幅大きくて違和感なしって感じ。
最後に
QQプロットをとりあえず作成するなら一瞬だけど、色々確認するとまだ勉強不足で時間かかりますな。