UNIX timestamp

まとめ

  • UNIX timestampのtimezoneはなにも指定しないとマシン設定
  • UNIXはミリ秒まで含んだ13桁のときがあるので注意。1000で割って10桁にするとdatetimeへ変換可能。

本文

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) in

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超えと大きいのでapplynp.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)

条件付き分布

標準正規分布・正規分布・多変量正規分布 | AGIRobots

多変量正規分布確率密度関数導出。数式は追えるけど、解釈はそれほど。

因子分析まとめ

自分用メモです。間違えていることあります。

発展因子分析 - YouTube

因子が解釈しやすいとは?スケールなし、0から正の範囲だと◎。それを数学的に表現すると2乗和の分散最大化。これがバリマトリックス回転(直交)。さらに直交を外して解釈しやすくしたいのでプロマックス回転。

重回帰分析参考サイト

自分用のメモです。言葉足らずですし、間違えていることもあります。

重回帰モデルの理論と実装 -なぜ正則化が必要か- | Deploy on Friday

正則化ってなにをしてるかについて。正規方程式で^tXXが正則でない(逆行列もたない)ときになんらかで調整したい。例えばリッジは誤差二乗和にλ| \boldsymbol{\beta} |^2を足すことで、\boldsymbol{\beta}微分したりして進めたときの正規方程式の左辺が ^tXX + {\lambda}Iみたいな感じになる。

幾何学でみる重回帰の性質 | 有意に無意味な話

幾何的にみた重回帰の解釈。 \boldsymbol{y}Xを作用させたら勝手にL(X)に射影してくれて、射影した先\hat{\boldsymbol{y}} -  \boldsymbol{y}L(X)は直交してくれるものだと思ってた。あまりまだ正しく理解できていない気がする。

正規方程式(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", # トキュメントにはあるのに使えなかった
);

f:id:iiiiikamirin:20220219141959p:plain

例の通り、各資産についてやってみた。自己相関がそんなに強いものはあまりない印象。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")

f:id:iiiiikamirin:20220219143639p:plain

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()

f:id:iiiiikamirin:20220219153003p:plain TPX版 f:id:iiiiikamirin:20220219153026p:plain

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")

f:id:iiiiikamirin:20220219153144p:plain

まとめ

株価等のデータは裾が厚いという特徴あり+相関は外れ値の影響受けやすい=全データを対象に自己相関について判断すると米株は逆相関が強いみたいなミスリードが生じるので注意が必要ですな。

特に株価が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使うのが一番簡単そう。なので、使ったうえで、パラメータの挙動諸々を確認してみる。

使用データ

ドル円の対数リターンを使います。 f:id:iiiiikamirin:20220209091712p:plain

実際のコード

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()

f:id:iiiiikamirin:20220209101149p:plain

stats.modelのドキュメントを見ながら、パラメータなどについて確認してみます。主だったパラメータは以下の通りだと思います。

  • dist:比較する分布(正規分布など)
  • distargs:比較する分布作成のためのパラメータ
  • fit:与えられた分布を平均や標準偏差で正規化するか

不変分散を基に母集団について検定するので、設定する分布は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()

f:id:iiiiikamirin:20220219132423p:plain

最後に作成したプロットを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)

f:id:iiiiikamirin:20220219131648p:plain

こうみると、上も(正規分布に属しているなら)もっと小さな値で上位数パーセントの数値が出てくるところを、より大きな値で上位数パーセントが出てしまっている=右側の裾があつい、と解釈。下側も同様に解釈して左右の裾が正規分布と比較してあつくなっていることをわかりやすく表せました。

ドル円以外でもQ-Qプロットを作成してみる

ドル円以外はどうなのでしょう??TPXは下落の方が裾が大きくなっているみたいだけど、SPXはそんなことなさそう。(下落時の方が値幅は大きくなりそうなので、SPXの結果は意外)。為替はドル円は左右でそんなに変わらないけどEURUSDは上昇の方が大きいのが意外。リスクオフではドル買いでは。債券は米債欧州債ともに売られる方が値幅大きくて違和感なしって感じ。 f:id:iiiiikamirin:20220209111755p:plain f:id:iiiiikamirin:20220209111821p:plain

最後に

QQプロットをとりあえず作成するなら一瞬だけど、色々確認するとまだ勉強不足で時間かかりますな。