指数分布による為替変動のシミュレーション・その2

こんばんは、皆様、三頌亭です。さて小ネタ集の続きですｗ。というか補足でございます。為替レートの変動（1分足）6カ月分をストレージしていきなり指数分布に当てはめました。こういった時系列変動がどのような分布を取るかについては様々な研究がありまして、やれ正規分布だとかべき乗分布、レヴィ分布とうるさいですｗ。ただ生データを見れば明らかなように正規分布などよりはウェストの細い、裾広がりの分布（「魔女の帽子」形）になっていることがわかると思います。つぎに生データの出現頻度をlogスケールにしたものと正規分布＆レヴィ分布と指数分布の比較のグラフをお示しいたします。下の2つのグラフは下記のところから拝借いたしました（株式です）。（Eugene Stanleyあたりの論文からだと思うのですがファイルが探せないです）

Statistical Properties of Stock Prices (Modeling the Dynamics of Stock Markets)

https://towardsdatascience.com/statistical-properties-of-stock-prices-15145be752a2

今回は大体18万個の生データを使いましたがとにもかくにも個人でこんなことができるようになったことに隔世の感ひとしおでございますｗ。すそ野がまだ安定してないのはデータ数の加減（出現頻度）でしょう。一般の短期のトレーダーさんにはあまり必要のない領域の出現頻度ではないかと三頌亭は思いますが・・・。

指数分布01

指数分布02

指数分布03

指数分布による為替変動のシミュレーション

こんばんは、皆様、三頌亭です。これも小ネタ集のひとつです。今回はVBAを使わなくてもできるでしょうｗ。といいますか1種のフィールドワークのようなものです。まず表題通り、為替の分足データをゲットしてきます。分足データはいろいろあるのですが下記のところで拝借いたしました。通貨ペアはドル円を選びました。

https://www.forexite.com/free_forex_quotes/forex_history.html

ここのデータで1分間の為替レートの変動を6カ月分ストレージします。変動は終値の差分を取りました。ここのところ対数正規だとかいろいろあるのですが、為替レートの場合株式などよりヴォラティリティがだいぶ低いので差分でもまあいいでしょうｗ。

まずこのデータをもとに度数分布を調べて、これを指数分布（ラプラス分布）に当てはめます。これをもとにエクセルの一様乱数から指数分布の乱数を発生させます。指数分布の乱数については下記が分かりやすいです。

http://godfoot.world.coocan.jp/Exp-Rand.htm

1分間ドル円の変動(λ=0.7,pips) = -(1/λ)*ln( RAND( ) )*(-1)^(RANDBETWEEN(1,10))

(正負両側にするため再度±1を乱数で発生させ乗算。単位をpipsからyenにするためには0.01を掛ける）

この指数分布の乱数を発生させて1日分の為替レートの変動をシミュレーションいたします。昔のデータなので111円をスタートにしてあります。発生した変動分を順次加算していくと1日の為替レートの変動をシミュレーションできることになります。

グラフ1～6番目はそのときの変動の1例です。テクニカル指標でよく使われるボリンジャーバンド（50分）を一緒に表示してあります。どうですか？なかなかそれらしいですねｗ。実は乱数を使っても、トレンドフォローにいいパターンや雇用統計のような急なパターンの変動が生まれるのかどうかを知りたかったのでこんなものを作ってみました。

因みに指数分布の性質に「無記憶性」というものがあります。わかりやすく言うとどの点でも次の1分後の変動の出現確率は同じで過去のデータには左右されないということです。これではトレーダーさんは「no hope」ですねｗ～困りましたｗ

頻度分布(USD/JPY)

指数分布当てはめ

シミュレーション例01

シミュレーション例02

シミュレーション例03

シミュレーション例04

シミュレーション例05

シミュレーション例06

連続データ重ね合わせ用・関数

こんばんは、皆様、三頌亭です。今日は小ネタ集の一つですｗ。以下の関数は縦に並んだデータを一定量シフトさせて足し合わせていく関数です。何に使うのかよくわからないという方が大半でしょうが、必要な人には自分でいうのもなんですが「便利」ですｗ。

Rem 連続データ重ね合わせ　配列として出力（配列数式使用）

Rem xtstp：繰り返しセルエリア（同一カラムで指定）

Rem deltax：シフトさせるセル数

Rem rep：シフトの繰り返し数

Rem addsh：シフトをスタートさせるセル位置（0が入力される）

Public Function multidoseadd(xtstp As Variant, deltax As Double, rep As Double, addsh As Double) As Variant

   Dim input1(8000, 1) As Variant

   Dim cnt As Double

   Dim result(8000, 1) As Variant

   cnt = 0

   deltax = Int(deltax)

   addsh = Int(Abs(addsh))

   Rem 指定範囲から順番にデータ用セルを取り出す（ゼロ時間から）

   For Each myCell In xtstp

      input1(cnt, 0) = myCell.Value

      cnt = cnt + 1

   Next

   Rem 重ね合わせ(deltax=シフト数）

    For n = 1 To rep

    For j = addsh To cnt - 1 + addsh

    result((n - 1) * deltax + j, 0) = result((n - 1) * deltax + j, 0) + input1(j - addsh, 0)

    Next j

    Next n

    multidoseadd = result

 End Function

注意書きです。

「Microsoft 365用の2018年9月の更新プログラムから、複数の結果を返すことができるすべての数式は、自動的に下方向に、または隣接するセルにスピルされます。 この動作の変更には、いくつかの新しい 動的な配列関数も伴います。 既存の関数または動的な配列関数のどちらを使用している場合でも、1つのセルに入力するだけで、 enterキーを押すことで確定する必要があります。 以前の従来の配列数式では、最初に出力範囲全体を選択してから、 Ctrl + Shift + Enter を押して数式を確認する必要があります。 一般的に、 CSE 数式と呼ばれています。」

連続データ重ね合わせ用・関数・使用例

数値的逆ラプラス変換（誤差評価にかえて・その2：Bessel関数編）

こんばんは、皆様、三頌亭です。前回はあまり得意でない矩形波を逆変換した例をお見せいたしましたので、今回はベッセル関数を例にとって細野法で数値的逆ラプラス変換行ってエクセルの関数値と比較してみましょう。ところで最近ではフリーの数式処理システム（maximaなど）もありますので、それらで逆ラプラス変換をやってみて時間次元の解析解が簡単にに求まるならば数値計算にはそれを使えばよろしいです。逆変換が簡単な初等関数の組み合わせで表現できない場合、近似値を求める手段として意味があると三頌亭は考えます。

さて、ベッセル関数計算用のコードは以下の通りです。

Rem 数値的逆ラプラス変換filt22(t) Hosonoのアルゴリズム

Rem atstp：分母　btstp：分子---ラプラス関数値の算出に必要なパラメータをセル参照

Rem 作例 F(s) = (sqrt(s^2+1)-s)^n/sqrt(s^2+1) 第一種ベッセル関数Jn(t) :n=a01(0)

Public Function filt22(time As Double, atstp As Variant, btstp As Variant, nts As Integer, p As Double) As Double

Dim ss(10000) As Complex

Dim x1L(10000) As Complex

Dim res(10000) As Double

Dim aa As Double

Dim n As Double

Dim cnt As Double

Rem パラメータ関連・複素数変数

Dim resigma As Double

Dim ared(100, 1) As Double

Dim bred(100, 1) As Double

Dim pi As Double

Dim a01(20) As Complex

Dim b01(20) As Complex

Dim ht(20) As Complex

Dim gt(20) As Complex

Dim ht01 As Complex

Dim ht02 As Complex

Dim ht03 As Complex

Dim ht04 As Complex

Dim gt01 As Complex

Dim gt02 As Complex

Dim gt03 As Complex

Dim gt04 As Complex

Dim ap(1000) As Double

Dim cp(1000) As Double

Rem 各種定数

pi = 3.14159265358979

Rem オイラー変換係数

cp(0) = 1

ap(p + 1) = 0

For q = 1 To p

    cp(q) = cp(q - 1) * (p + 1 - q) / q

Next q

For q = p To 0 Step -1

    ap(q) = ap(q + 1) + cp(q)

Next q

Rem パラメータの取り込み・複素数変換

Rem 指定範囲から順番にデータ用セル値を取り出す（次数の高い順a0-a8,b0-b8）

For Each myCell In atstp

      ared(cnt, 0) = myCell.Value

      a01(cnt) = ToComplex(ared(cnt, 0), 0)

      cnt = cnt + 1

Next

cnt = 0

For Each myCell In btstp

      bred(cnt, 0) = myCell.Value

      b01(cnt) = ToComplex(bred(cnt, 0), 0)

      cnt = cnt + 1

Next

Rem 高速ラプラス変換（細野法）

If time > 0 Then

aa = 10

For n = 1 To nts + p

    ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)

    Rem  ラプラス変換関数記述　Cadd Csub Cmul Cdiv ToComplex s = ss(n): 

    gt01 = csqr(Cadd(Cmul(ss(n), ss(n)), ToComplex(1, 0)))

    ht01 = Cpow(Csub(gt01, ss(n)), a01(0))

    x1L(n) = Cdiv(ht01, gt01)

    Rem  ラプラス変換関数記述　Cadd Csub Cmul Cdiv ToComplex s = ss(n)

res(n) = Imz(x1L(n)) * (-1) ^ n

Next

Rem nts項分加算

resigma = 0

For n = 1 To nts

    resigma = resigma + res(n)

Next

Rem オイラー変換-p項分追加

For i = 1 To p

    resigma = resigma + res(nts + i) * ap(i) / 2 ^ (p)

Next

Rem 変換値出力

   filt22 = resigma * Exp(aa) / time

Else

Rem t=0の関数値：初期値の出力

   filt22 = 0

End If

End Function

Rem 複素数計算用vba-code

Rem https://kms130.hatenablog.com/entry/2020/12/03/232117

エクセルの内蔵の関数BESSELJでの計算と比較して、オイラー変換10項、未変換30項で相対誤差は10＾(-7)から10^(-8)程度です。０次から8次ぐらいまでやってみましたが非常に性能のいい近似式です。ｔ-＞∞できれいに収束するような関数に対しては素晴らしく近似精度がよろしいです。因みに↑のコードのaa は 8にしてあってだいたい有効桁7,8桁ということでほぼ注文通り収まってます（グラフでは重なってわかりません）。ラプラス変換した関数の記述部分は以下の3行のみです。

Rem F(s) = (sqrt(s^2+1)-s)^n/sqrt(s^2+1)： 第一種ベッセル関数Jn(t)

gt01 = csqr(Cadd(Cmul(ss(n), ss(n)), ToComplex(1, 0)))

ht01 = Cpow(Csub(gt01, ss(n)), a01(0))

x1L(n) = Cdiv(ht01, gt01)

ベッセル関数　オイラー変換10項、未変換30項

数値的逆ラプラス変換（誤差評価にかえて）

こんばんは、皆様、三頌亭です。細野法でオイラー変換の項数と未変換の項数を自由に設定できるようにしたVBAのコードです。これでおそらくラプラス変換した関数部分の記述を変えるだけでオッケーでしょうＷ。これらの数値的逆ラプラス変換法のアルゴリズムがあまり得意ではない方形波パルスの逆変換を試してみます。オイラー変換の項数と未変換の項数は20,20のものと80,400のもののグラフをお示しいたします。

Rem 数値的逆ラプラス変換filt10(t) Hosonoのアルゴリズム----------------

Rem atstp：分母 btstp：分子-ラプラス関数値の算出に必要なパラメータをセル参照

Rem 作例 F(s) = exp(-a*s)/s*(1-exp(-T*s)) 方形波パルス

Rem パラメータ入力：a01(0)=a,b01(0)=T

Public Function filt10(time As Double, atstp As Variant, btstp As Variant, nts As Integer, p As Double) As Double

Dim ss(5000) As Complex

Dim x1L(5000) As Complex

Dim res(5000) As Double

Dim aa As Double

Dim n As Integer

Dim cnt As Double

Rem パラメータ関連・複素数変数

Dim resigma As Double

Dim ared(100, 1) As Double

Dim bred(100, 1) As Double

Dim pi As Double

Dim a01(20) As Complex

Dim b01(20) As Complex

Dim ht(20) As Complex

Dim gt(20) As Complex

Dim ht01 As Complex

Dim ht02 As Complex

Dim ht03 As Complex

Dim ht04 As Complex

Dim gt01 As Complex

Dim gt02 As Complex

Dim gt03 As Complex

Dim gt04 As Complex

Dim ap(200) As Double

Dim cp(200) As Double

Rem 各種定数

pi = 3.14159265358979

Rem オイラー変換係数

cp(0) = 1

ap(p + 1) = 0

For q = 1 To p

    cp(q) = cp(q - 1) * (p + 1 - q) / q

Next q

For q = p To 0 Step -1

    ap(q) = ap(q + 1) + cp(q)

Next q

Rem パラメータの取り込み・複素数変換

Rem 指定範囲から順番にデータ用セル値を取り出す（次数の高い順a0-a8,b0-b8）

For Each myCell In atstp

      ared(cnt, 0) = myCell.Value

      a01(cnt) = ToComplex(ared(cnt, 0), 0)

      cnt = cnt + 1

Next

cnt = 0

For Each myCell In btstp

      bred(cnt, 0) = myCell.Value

      b01(cnt) = ToComplex(bred(cnt, 0), 0)

      cnt = cnt + 1

Next

Rem 高速ラプラス変換（細野法）

If time > 0 Then

aa = 8

For n = 1 To nts + p

    ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)

    Rem  ラプラス変換関数記述　Cadd Csub Cmul Cdiv ToComplex F(s : ss(n))

    ht01 = Cdiv(ToComplex(1,0),ss(n))

    gt02 = Csub(ToComplex(0,0),a01(0))

    gt02 = Cmul(ss(n),gt02)

    gt01 = Csub(ToComplex(1,0),Cexp(gt02))

    gt03 = Csub(ToComplex(0, 0), b01(0))

    gt03 = Cmul(ss(n), gt03)

    x1L(n) = Cmul(Cmul(ht01, gt01), Cexp(gt03))

    Rem  ラプラス変換関数記述　Cadd Csub Cmul Cdiv ToComplex  F(s : ss(n))

res(n) = Imz(x1L(n)) * (-1) ^ n

Next

Rem nts項分加算

resigma = 0

For n = 1 To nts

    resigma = resigma + res(n)

Next

Rem オイラー変換-p項分追加

For i = 1 To p

    resigma = resigma + res(nts + i) * ap(i) / 2 ^ (p)

Next

Rem 変換値出力

   filt10 = resigma * Exp(aa) / time

Else

Rem t=0の関数値：初期値の出力

   filt10 = 0

End If

End Function

どうでしょう。まあ相変わらず無駄の多いみっともないコードで申し訳ないです（＾＾；）。ギブズ現象がはっきりわかります。時間を大きくしていくともっとはっきりしてきます。私の場合、拡散方程式でしたので急激な変動があまりありません。なので、あんまりわかりませんでした.。ルンゲクッタ法で近似解を求めても解はほとんど同じでした。因みによく論文なんかで使われているベッセル関数なんかは大体大丈夫です。ご参考までに・・・。

追記：時間大きくすると誤差が大きくなるという意味がわかりづらいと思いましたので、グラフの例を二つ追加いたしました。矩形波の連続パルスです。

作例 F(s) = (1/s)*(1/(1+exp(-T*s))) 方形波パルス（周期波）で関数部分のコードは

ht01 = Cdiv(ToComplex(1, 0), ss(n))
gt01 = Cadd(ToComplex(1, 0), Cexp(Cmul(ss(n), Csub(ToComplex(0, 0), a01(0)))))
gt01 = Cdiv(ToComplex(1, 0), gt01)

x1L(n) = Cmul(ht01, gt01)

です。

オイラー変換20項、未変換20項

オイラー変換80項、未変換400項

オイラー変換20項、未変換40項

オイラー変換80項、未変換1000項

数値的逆ラプラス変換（Gaver-Stehfest 法について）

こんばんは、皆様、三頌亭です。このまえから解説してます逆ラプラス変換方法ですが、ご参考までに細野法ではない、Gaver-Stehfest 法について簡単にご説明いたします。まあ大体、利用者の立場としては「能書きはいいから早くプログラムをよこせ！」というのがよくあるお話ですねｗｗ。

ところでGaver-Stehfest 法についてですが複素数を必要とせず、実数のみで計算が完了いたします。ラプラス変換済みの関数F(s)にそのまま実数s=Ln(2)*n/time(n=1～20程度まで)を入力してF(s)の関数値を計算いたします。それぞれの値に定まった係数をかけて加算していきます。総和にLn(2)/timeを掛けて、timeにおける時間次元の関数値f(time)とします。因みにかける係数は下記のとおりです。一番下の画像に一般式をお示します。

vv(1) = -5.51146384479718E-06

vv(2) = 0.152386463844797

vv(3) = -117.465476190476

vv(4) = 17342.4493386243

vv(5) = -922806.928902116

vv(6) = 23774087.7871031

vv(7) = -349421166.19537

vv(8) = 3241369852.23187

vv(9) = -20276948307.2378

vv(10) = 89464829823.7973

vv(11) = -287020921147.103

vv(12) = 682992010281.513

vv(13) = -1219082330054.38

vv(14) = 1637573800842.02

vv(15) = -1647177486836.12

vv(16) = 1221924554444.23

vv(17) = -648806558817.535

vv(18) = 233316653213.707

vv(19) = -50913800705.4676

vv(20) = 5091380070.54676

値をご覧になればお分かりになると思いますが、振動してどんどん大きくなりますｗ。なのでｎを大きくとると相対誤差が非常に悪くなってしまいます。一応、一般のｎについて係数を計算するvbaのコードを紹介いたします（金融工学の論文からのコピーです。貼り付ければそのままで動作します）。

 https://www.researchgate.net/publication/228310935_Option_Pricing_with_Jumps

 一番最後の部分がそうです。

Rem StehfestCoefficients returns coefficients in Stehfest inversion formula
Function StehfestCoefficients(Nterms As Integer) As Variant
Dim i As Integer, j As Integer, k As Integer, nh As Integer, sn As Integer, nn As Integer
nh = Nterms / 2 'Must be even
Nterms = 2 * nh
ReDim Coeffs(0 To Nterms) As Double, g(0 To Nterms) As Double, h(0 To nh) As Double
g(0) = 1
For i = 1 To Nterms
g(i) = g(i - 1) * i
Next i
h(1) = 2# / g(nh - 1)
For i = 2 To nh
h(i) = Exp(nh * Log(i)) * g(2 * i) / (g(nh - i) * g(i) * g(i - 1))
Next i
If nh Mod 2 = 0 Then sn = 1 Else sn = -1
For i = 1 To Nterms
Coeffs(i) = 0
If i < nh Then nn = i Else nn = nh
j = Int( (i + 1) / 2)
For k = j To nn
Coeffs(i) = Coeffs(i) + h(k) / (g(i - k) * g(2 * k - i))
Next k
sn = -sn
Coeffs(i) = Coeffs(i) * sn
Next i
StehfestCoefficients = Coeffs
End Function

とりあえず、正確に動作するようです。一般に発散していく関数については数値的逆ラプラス変換は不正確な値になりやすいです。

参考書（＾＾；）

 ALAN M. COHEN 「NUMERICAL METHODS FOR LAPLACE TRANSFORM INVERSION」
http://fulviofrisone.com/attachments/article/462/Springer%20-%20Numerical%20Methods%20for%20Laplace%20Transform%20Inversion%20-%202007.pdf

追記：因みに項数nはどのぐらいがいいのかという問題があるのですが、上記の論文では10-20の偶数がいいそうです。いくつか使ってる論文当たってみたのですが14か16とかぐらいが多いようです。わたしはめんどくさいので20個に固定して先計算した値を上記のようにして貼り付けて使ってましたｗ。

Gaver-Stehfest 法

Stehfest table CALCULATION-FORMULA

数値的逆ラプラス変換(vba-code解説)

こんばんは、皆様、三頌亭です。今日は数値的な逆ラプラス変換をエクセルに関数として組み込むプログラムを紹介いたします。もともとこのプログラムは実験データを拡散方程式の解に当てはめるために作成したものです。もうだいぶ前に作成したものですが、最近になっても取っつきやすいエクセル上で使用できるものがないようなので公開いたします。とはいうもののそんな大げさな、かっこいいものではありません（＾＾；）。しかし実戦使用でだいぶ使ってきましたので、愛着のあるプログラムですｗ。

ところで複素数計算のコードは下記のhpよりお借りして、よく使う双曲線関数と平方根の関数を追加して作成してあります。

Fallen Physicist, Rising Engineer

https://sci.tea-nifty.com/blog/2008/10/excelvba-3b05.html?optimized=0

まずこのパッケージをエクセルから開発＞visual basic＞挿入＞標準モジュールで貼り付けてください。そののちに逆変換のfilt06を追加で貼り付けてxlmsで保存して出来上がりです。

filt06(time As Double, atstp As Variant, btstp As Variant, nts As Integer) となっております。関数計算に必要なパラメータ・セットはatstpとbtstpでセル範囲を指定します。作例では有理関数の分母分子のsの8次式までの係数です。Ntsはオイラーの収束加速をもちない項数を入力します。20-40ぐらいでしょうか？。オイラーの収束加速は10項で固定してありますが、下記の係数のプログラムを組み込むことで大きくすることは可能です。

Public function eulercoeff(p as double, eff as double) as double

dim ap(100) as double

dim cp(100) as double

cp(0)=1

ap(p+1)=0

for q=1 to p

    cp(q)=cp(q-1)*(p+1-q)/q

next q

for q=p to 0 step -1

    ap(q)=ap(q+1)+cp(q)

next q

eulercoeff = ap(eff)

end function

数値的逆ラプラス変換のアルゴリズムはたくさんのものがあります。Stehfestのものをはじめとしていくつかのものを作成してみましたが、細野のものが最もロバストだったので現在もこれを使用しています。もちろん、時間と空間のメッシュを作成して実関数のみを使用する近似解法（ルンゲクッタ、クランク-ニコルソン）もあるのですが実験データの当てはめ（非線形最小二乗法）との相性があまりよくないのでこの方法に落ち着きました。ご参考までに追記しておきます。

参考文献

BASICによる高速ラプラス変換（細野敏夫：共立出版：1984）

filt001

スクリーンショット

filt002