水曜日のドル円は上昇する
- 曜日によってドル円の上昇に差はあるのかを調べてみる。
- いつものようにプロットまで一気に。
import pandas as pd from pandas.plotting import autocorrelation_plot import numpy as np import scipy from scipy import stats as st from matplotlib import pylab as plt import seaborn as sns sns.set() import statsmodels.api as sm from statsmodels.stats.multicomp import pairwise_tukeyhsd import heapq as hp import math import datetime import time fx_data = pd.read_csv("/Users/ユーザー名/Documents/FX/foreign_exchange_historical_data/USDJPY/USDJPY_DAY.csv")#データの読み込み #pandasの計算は遅いので、全てnumpy配列に変換 date = np.array(fx_data["date"]) opening = np.array(fx_data["opening"])#始値 high = np.array(fx_data["high"])#高値 low = np.array(fx_data["low"])#低値 closing = np.array(fx_data["closing"])#終値 dict = {"Monday":0, "Tuesday":0, "Wednesday":0, "Thursday":0, "Friday":0}#値上がりする日の集計 dict_all = {"Monday":0, "Tuesday":0, "Wednesday":0, "Thursday":0, "Friday":0}#曜日ごとの日数 for i in range(len(date)):#dictもdict_allも一度に計算 temp = datetime.datetime.strptime(date[i], "%Y/%m/%d") temp = temp.strftime("%A") dict_all[temp] += 1 if closing[i] >= opening[i]: dict[temp] += 1 ratio_of_increase = [a/b for a, b in zip(dict.values(), dict_all.values())] #割合に変換 # %%まとめて選択 plt.bar(dict.keys(), ratio_of_increase, color = "skyblue", alpha = 0.7) plt.title("usdjpy") plt.xlabel("week") plt.ylabel("increase") plt.savefig("/Users/ユーザー名/Desktop/week.png", format = "png", dpi = 300)
- 図示はこちら。
- 水曜日が若干上昇しやすい印象もあるが・・・。
- χ二乗検定にかけてみる。
dict_exp = [c//2 for c in dict_all.values()]#χ二乗検定の期待値は0.5とする。 dict_copy = [d for d in dict.values()] stats.chisquare(dict_copy, f_exp = dict_exp)#χ二乗検定 dict_nega = [f-e for e, f in zip(dict.values(), dict_all.values())] st.chi2_contingency(np.array([dict_copy, dict_nega]))#tableでも同様に
結果は以下。
Power_divergenceResult(statistic=1.2519046672629104, pvalue=0.8694811382649527) (2.546521721041146, 0.63632517975807, 4, array([[347.39930955, 347.89988493, 348.4004603 , 348.4004603 , 347.89988493], [346.60069045, 347.10011507, 347.5995397 , 347.5995397 , 347.10011507]]))
- 有意差はないようだ。
- 解析期間(2007年4月2日から2020年8月15日）
dplyrの使い方（コードの簡略化）
- この記事でのforの繰り返しをdplyrの関数で簡略化してみる。
- 下記の記事を参考にした。
- Rのソースコードはこちら。
library(ggplot2) library(dplyr) library(RcppRoll) d <- read.csv("~/Documents/FX/USDJPY/USDJPYD.csv") N <- nrow(d) #データ数 D <- 245 #日数の最大値の設定 z1 <- array(rep(0, D)) #空の配列を用意、後で割合の集計に利用 z2 <- array(rep(0, D)) #空の配列を用意、後で割合の集計に利用 for (i in 1:D){ d <- d %>% #ローリング関数で区間中の最大値・最小値を計算 dplyr::mutate(rmax = roll_max(d$high, i, weights = NULL, align = "left", fill = NA, na.rm = FALSE)) %>% dplyr::mutate(rmin = roll_min(d$low, i, weights = NULL, align = "left", fill = NA, na.rm = FALSE)) %>% #条件設定 dplyr::mutate(ask = if_else(high < rmax, 1, 0)) %>% dplyr::mutate(bid = if_else(low > rmin, 1, 0)) %>% #累積関数 dplyr::mutate(sum_ask = cumsum(ask)) %>% dplyr::mutate(sum_bid = cumsum(bid)) z1[i] <- d$sum_ask[N-i]/(N-i) z2[i] <- d$sum_bid[N-i]/(N-i) } z <- data.frame(day = 1:D, z1, z2) #データフレームの作成 ggplot() + theme_set(theme_bw(base_size = 14)) g <- ggplot(z, aes(x = z$day)) g <- g + geom_line(aes(y = z$z1, colour = "ask"), size = 1) g <- g + geom_line(aes(y = z$z2, colour = "bid"), size = 1) g <- g + xlab("day") + ylab("probability") g <- g + scale_x_continuous(breaks = seq(0, D, by = 30)) g <- g + labs(colour = "ask/bid") g <- g + scale_y_continuous(breaks = seq(0.4, 1, by = 0.1), limits = c(0.4, 1)) plot(g)
- プロットは省略。結局for文を使わずにコーディングすることは難しいか・・・。
- ポイントは累積関数を用いているところ。なぜか
z1[i] <- sum(d$ask == 1)/(N-i)
とするとエラーが出て進まなかった。
- 二重for文よりは遥かに早い。
変動値の従う分布(8)
data { #データYの宣言 int N; vector[N] Y; } parameters { #サンプリングしたいパラメータθの宣言 real<lower = 0> alpha; real<lower = 0> theta; } model { #尤度p(Y|θ)の記述 #事前分布p(θ)の記述 #指数分布を仮定 Y ~ gamma(alpha, theta); }
- R kickはこちら。
d <- read.csv("~/Documents/FX/USDJPY/USDJPYD.csv") co <- d$closing-d$opening ho <-d$high-d$opening ol <- d$low-d$opening oo <- diff(d$opening, lag = 1) d <- data.frame(N = nrow(d), co, ho, ol) N <- nrow(d) data <- list(N = nrow(d), Y = d$ho) fit <- stan(file = "~/Documents/FX/USDJPY/Rdata/distribution_model2.stan", data = data, seed = 1234) ms <- rstan::extract(fit) new_data <- data.frame(ho) alpha <- mean(ms$alpha) lambda <- 1/mean(ms$theta) myfun <- function(x){ 1/gamma(alpha)/(lambda)^(alpha)*x^(alpha-1)*exp(-x/lambda) } ggplot() + theme_set(theme_bw(base_size = 14)) g <- ggplot(new_data) g <- g + geom_histogram(aes(x = new_data$ho, y = ..density.., fill = "high-opening"), binwidth = 0.1) g <- g + labs(fill = "high-opening") g <- g + stat_function(fun = myfun, colour = c("#EE3B3B"), size = 1 ) g <- g + xlab("high-opening") plot(g)
- 結果と図示はこちら。なかなかいい感じではないだろうか。
nference for Stan model: distribution_model2. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff alpha 1.33 0.00 0.03 1.27 1.31 1.33 1.35 1.38 1248 lambda 3.05 0.00 0.08 2.89 2.99 3.05 3.11 3.21 1256 lp__ -434.06 0.02 0.95 -436.53 -434.48 -433.78 -433.36 -433.11 1587 Rhat alpha 1 lambda 1 lp__ 1
- （始値-安値）の方は結果のみ。
Inference for Stan model: distribution_model2. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff alpha 1.22 0.00 0.03 1.17 1.20 1.22 1.24 1.28 1078 theta 2.48 0.00 0.07 2.35 2.43 2.48 2.53 2.63 1090 lp__ -847.26 0.03 1.05 -850.10 -847.67 -846.92 -846.51 -846.23 1242 Rhat alpha 1 theta 1 lp__ 1
変動値の従う分布(7)
data { //データYの宣言 int N; vector[N] Y; } parameters { //サンプリングしたいパラメータθの宣言 real df; real mu; real sigma; } model { //尤度p(Y|θ)の記述 //事前分布p(θ)の記述 Y ~ student_t(df, mu, sigma); }
- これをRでキック。プロットまで一気に行う。t分布の密度関数は自分でセコセコ打ち込まなくてはいけないのが悲しい*1。
d <- read.csv("~/Documents/FX/USDJPY/USDJPYD.csv") co <- d$closing-d$opening ho <-d$high-d$opening ol <- d$low-d$opening oo <- diff(d$opening, lag = 1) d <- data.frame(N = nrow(d), co, ho, ol) library(rstan) N <- nrow(d) data <- list(N = nrow(d), Y = d$co) fit <- stan(file = "~/Documents/FX/USDJPY/Rdata/distribution_model.stan", data = data, seed = 1234) ms <- rstan::extract(fit) new_data <- data.frame(co) v <- mean(ms$df) mu <- mean(ms$mu) sigma <- mean(ms$sigma) myfun <- function(x){ gamma((v+1)/2)*(1+(x-mu)^2/(v*(sigma)^2))^(-(v+1)/2)/(sigma*gamma(v/2))/sqrt(v*(pi)) } ggplot() + theme_set(theme_bw(base_size = 14)) g <- ggplot(new_data) g <- g + geom_histogram(aes(x = new_data$co, y = ..density.., fill = "closing-opening"), binwidth = 0.1) g <- g + labs(fill = "closing-opening") g <- g + stat_function(fun = myfun, colour = c("#EE3B3B"), size = 1 ) g <- g + xlab("high-opening") plot(g)
- 結果は以下。当てはまりは悪くはないが、尖度が今ひとつといったところ。
- 拡大してみる。
- より拡大するとやはり単純にt分布に従うとは言えないと考えられる。
Inference for Stan model: distribution_model. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff df 3.65 0.01 0.27 3.16 3.46 3.64 3.83 4.23 2259 mu 0.01 0.00 0.01 -0.01 0.00 0.01 0.01 0.03 2651 sigma 0.47 0.00 0.01 0.44 0.46 0.47 0.47 0.49 2310 lp__ -1125.89 0.03 1.22 -1128.99 -1126.46 -1125.57 -1124.98 -1124.50 1878 Rhat df 1 mu 1 sigma 1 lp__ 1
*1:なにか良い方法があれば追記
変動値の従う分布(6)
- ドル円の変動値が従う分布を推定してみる。t分布や指数分布のパラメータは固定するものとする。
- まずは（高値-始値）および（始値-安値）の従う指数分布のパラメータ
から推定してみる。
- 推定にはrstanを用いる。
data { #データYの宣言 int N; vector[N] Y; } parameters { #サンプリングしたいパラメータθの宣言 real beta; } model { #尤度p(Y|θ)の記述 #事前分布p(θ)の記述 #指数分布を仮定 Y ~ exponential(beta); }
- 読み込みは次の様になる。
d <- read.csv("~/Documents/FX/USDJPY/USDJPYD.csv") co <- d$closing-d$opening ho <-d$high-d$opening ol <- d$low-d$opening oo <- diff(d$opening, lag = 1) d <- data.frame(N = nrow(d), co, ho, ol) library(rstan) N <- nrow(d) data <- list(N = nrow(d), Y = d$ho) fit <- stan(file = "~/Documents/FX/USDJPY/Rdata/distribution_model.stan", data = data, seed = 1234) ms <- rstan::extract(fit)
- 結果は以下。
mean se_mean sd 2.5% 25% 50% 75% 97.5% beta 2.30 0.00 0.04 2.21 2.27 2.30 2.33 2.38 lp__ -503.01 0.02 0.74 -505.16 -503.19 -502.73 -502.53 -502.48 n_eff Rhat beta 1468 1 lp__ 1886 1
- きちんと収束している。図示は次。
myfun <- function(x){ mean(ms$beta)*exp(-mean(ms$beta)*x) } new_data <- data.frame(ho, ol) ggplot() + theme_set(theme_bw(base_size = 14)) g <- ggplot(new_data) g <- g + geom_histogram(aes(x = new_data$ho, y = ..density.., fill = "high-opening"), binwidth = 0.1) g <- g + labs(fill = "high_low") g <- g + stat_function(fun = myfun, colour = c("#EE2C2C"), size = 1 ) plot(g)
- ggplot2でヒストグラムを描写する時は、countとdensityの2種類のどちらかを指定できるが、ここではdensityを指定しておかないとうまく重ね書きができない。
- 当てはまりは良好そうである。（始値-安値）でも推定すると、
となる。
変動値の従う分布(5)
- 分布を調べるのに今度はghypパッケージを使用してみる。
- AICをもとに一般化双曲型分布、normal inverse Gaussian分布、Variance Gamma分布、歪んだt分布などへの当てはまりを調べることができる。
https://cran.r-project.org/web/packages/ghyp/ghyp.pdf
d <- read.csv("~/Documents/FX/USDJPY/USDJPYD.csv") co <- d$closing-d$opening ho <-d$high-d$opening ol <- d$low-d$opening oo <- diff(d$opening, lag = 1) stepAIC.ghyp(co)
- まずはstepAICで最も当てはまりの良い分布を探索。
$fit.table model symmetric lambda alpha.bar mu sigma 8 NIG TRUE -0.500000000 0.7295650 0.002865157 0.6549768 7 hyp TRUE 1.000000000 0.2404656 -0.002839829 0.6481153 9 VG TRUE 1.140484515 0.0000000 -0.006562022 0.6470659 6 ghyp TRUE 0.130331659 0.6745057 0.001558661 0.6522863 3 NIG FALSE -0.500000000 0.7326138 0.012276449 0.6546704 4 VG FALSE 1.133472700 0.0000000 -0.011106709 0.6473282 2 hyp FALSE 1.000000000 0.2356592 -0.005071144 0.6482965 1 ghyp FALSE 0.008506488 0.7012681 0.008752815 0.6525317 10 t TRUE -1.821723585 0.0000000 0.004649114 0.6809921 5 t FALSE -1.831107876 0.0000000 0.017310757 0.6796988 11 gauss TRUE NA Inf -0.001664452 0.6555834 gamma aic llh converged n.iter 8 0.000000000 5573.093 -2783.546 TRUE 122 7 0.000000000 5573.315 -2783.658 TRUE 146 9 0.000000000 5573.386 -2783.693 TRUE 78 6 0.000000000 5574.444 -2783.222 TRUE 287 3 -0.013992090 5574.642 -2783.321 TRUE 149 4 0.009427421 5574.985 -2783.492 TRUE 373 2 0.003389183 5575.284 -2783.642 TRUE 233 1 -0.010626758 5576.207 -2783.104 FALSE 502 10 0.000000000 5582.168 -2788.084 TRUE 86 5 -0.019432190 5583.372 -2787.686 TRUE 193 11 0.000000000 6003.187 -2999.593 TRUE 0
- NIGはnormal inverse Gaussian分布のこと。
Normal-inverse Gaussian distribution - Wikipedia
- 高値、安値でも同様にstepAICを行う。結果は省略するが、optim（尤度関数への当てはめ）が収束している分布の中ではやはりNIGが最も上位にあがる。
- （翌日の始値-当日の始値）ではVariance Gammaが最上位に。ただし、NIGも3番目と悪くない。
- データのまとめも簡単で、以下のコードで済む*1。
fit1 <- fit.NIGuv(co) summary(fit1) hist(fit1)
- 結果はこちら。
Asymmetric Normal Inverse Gaussian Distribution: Parameters: alpha.bar mu sigma gamma 0.73261378 0.01227645 0.65467042 -0.01399209 Call: fit.NIGuv(data = co) Optimization information: log-Likelihood: -2783.321 AIC: 5574.642 Fitted parameters: alpha.bar, mu, sigma, gamma; (Number: 4) Number of iterations: 149 Converged: TRUE
*1:例えば分布に対称性を課す場合は、symmetric = TRUE とすれば良い
