【Rで非ガウス統計４】間欠性 (Intermittency)と分散不均一性 (Heteroscedasticity)：多重スケール確率分布解析（その１）

私たちは日常の中で、「変化の少ない期間が続いたかと思うと、突然、大きな変動が立て続けに起こる」という経験をすることがあります。静かな天気が続いていたのに、急に激しい嵐に見舞われることもあれば、株式市場が落ち着いていると思っていた矢先に、ある日を境に価格が大きく揺れ動き始めることもあります。このような現象は、実は多くの分野で共通して観測されます。
　長年にわたり、物理学者たちを悩ませてきた問題の一つが乱流です。乱流とは、川の急流や煙の立ちのぼりに見られるような、不規則で複雑な流れの状態を指します。見た目には完全にランダムに見えるこの流れも、統計的な視点で眺めると、意外な秩序が浮かび上がってきます。
　もし乱流が空間的にも時間的にも均一であれば、速度の分布は平均値のまわりに広がる正規分布になるはずです。ところが実際の乱流では、大きなスケールで見たときにはガウス分布に近い形を示す一方で、観測スケールを小さくするにつれて、非常に大きな速度変動が起こる確率が高くなり、裾の厚い非ガウス分布が現れます。
　このような分布の変化が生じるのは、乱流中のエネルギーが一様に分布していないためです。穏やかな流れが続く領域の中に、細かな渦や急激な速度変化が局所的に集中する領域が点在しており、観測点がそうした領域を通過した瞬間に、極端な速度変動が捉えられます。その結果、乱流全体としては同じ強さに見えても、局所的には静かな状態と激しい状態が混在することになります。
　このように、乱流の速度分布が正規分布から外れるという事実は、激しい運動が間欠的に現れるという乱流固有の性質を反映しており、この性質は間欠性（intermittency）と呼ばれています。
　一方、金融市場の株価指数や為替レートを眺めていると、こちらにも乱流とよく似た光景が見られます。値動きが小さい静かな期間がしばらく続いたかと思うと、突然、価格が大きく上下する時期が訪れます。そして不思議なことに、その荒れた状態は一日で終わることは少なく、しばらく続く傾向があります。経済学や統計学の分野では、このように「変動の大きさそのものが時間とともに変化する」性質を、分散の不均一性（heteroscedasticity）と呼びます（下図参照）。

分散が一定な時系列（上）と分散の不均一な時系列（下）の例。Homoscedasticity は「等分散性」、Heteroscedasticity は「分散不均一性」を意味する。灰色の領域は局所標準偏差の１倍（内側）と２倍（外側）の幅を表す。

　金融の専門家は、この現象を「市場が不安定になった」「ボラティリティが高まった」と表現します。しかし、統計的な視点から眺めると、ここで起きていることは乱流の間欠性と非常によく似ています。すなわち、価格変動の分散が常に一定なのではなく、静かな時期と荒れた時期が共存しているのです。ここでもまた、平均や分散を一つ計算しただけでは、現象の本質を捉えることはできません。
　このように見てくると、乱流と金融市場は、表面的にはまったく異なる世界でありながら、「ゆらぎの強さが一様ではない」という共通の特徴をもっていることが分かります。物理学ではそれを「間欠性」と呼び、金融や統計の分野では「分散の不均一性」と呼ぶ――その違いは主に用語の違いにすぎません。名前は異なっていても、見ている現象の本質は同じです。強い変動はいつでも起こるわけではなく、特定の時間や状況に集中して現れ、その結果として、極端な出来事が無視できない頻度で生じます。
　この共通点に気づくと、自然と次の疑問が浮かび上がります。「では、こうした現象をどのように定量的に捉えればよいのだろうか」。平均値や分散だけでは不十分であるなら、別の視点が必要になります。その一つの答えが、「スケールを変えながら揺らぎの分布を見る」という考え方です。小さな時間幅で見たときと、大きな時間幅で見たときとで、揺らぎの姿はどのように変わるのか。その変化の仕方自体に、間欠性や分散の不均一性の情報が刻み込まれています。
　この発想から生まれたのが、多重スケールで確率分布を眺めるという解析の枠組みです。乱流の渦も、株式市場の価格変動も、そして心拍のゆらぎも、異なるスケールで見直すことで、これまで見えなかった秩序や構造を語り始めます。ということで、今回は「多重スケール確率分布解析」を紹介します。とはいえ、詳しく説明するには、長い説明が必要ですので、まずは「分散の不均一性」と「粗視化」についてイメージがつかめるように解説します。

1. 多重スケール確率分布解析
2. 独立同分布の非ガウス過程と間欠性過程の違い

1. 多重スケール確率分布解析

　時系列データの解析では、平均値や分散、あるいはパワースペクトルなどがよく用いられます。しかし、乱流、金融データ、心拍変動のような複雑な時系列では、それらの指標だけでは捉えきれない特徴が現れることがあります。その代表例が、非ガウス性と間欠性です。
　非ガウス性とは、確率分布が正規分布よりも裾の厚い形を示す性質を指します。一方、間欠性とは、大きな変動が局所的に集中して現れ、その合間に小さな変動が続くといった、「分散の非一様性」をもつ性質です。
　多重スケール確率分布解析は、このような性質を、「時系列を観察する解像度（粗視化スケール）を変えながら確率分布を見る」ことで定量化しようとする手法です。従来のように平均値や分散を一つ計算する解析とは異なり、粗視化スケールを変化させつつ、確率分布そのものの変形を追跡する点に大きな特徴があります。下図は「フルーツ」の画像を粗視化したもので、粗視化とは「解像度を低くする処理」です。
　以下では、多重スケール確率分布解析の流れを順に説明します。

粗視化のイメージ。画像の解像度を粗くするように、時系列を粗視化（局所平均化）して観察するのが多重スケール解析。

1-1. 粗視化（coarse-graining）

　多重スケール確率分布解析ではまず、平均値0、有限分散をもつ定常時系列 $\displaystyle{ \{x _ i\} }$ を考えます。この時系列に対し、粗視化スケールを $s$ とし、連続する $s$ 点の和をとります。

$\displaystyle{ \Delta_s z_i = \sum_{j=1}^{s} x_{i+j} }$

これは、積分時系列

$\displaystyle{ z_i = \sum_{k=1}^{i} x_{k} }$

の差分と等価であり「時間スケール $s$ で見た変動量」に対応します。あるいは、 $\displaystyle{ \Delta _ s z _ i }$ を $s$ で割ると局所平均になりますので、スケールを変えて局所平均のばらつきを観察すると考えても構いません（上図では画像を平均値で粗くしています）。なぜなら、我々の興味は分布の形にあるので、 $s$ で割っても割らなくても結果（分布の形）は変わりません。

1-2. 標準化と確率分布

　粗視化した時系列 $\displaystyle{ \Delta _ s z _ i }$ を平均0、分散1に標準化し、その確率密度関数（PDF）を推定します。このPDFの形の変化を、スケール $s$ を変えながら調べていきます。

1-3. 非ガウス分布の定量化

　ここでは、標準化された粗視化時系列 $\displaystyle{ \Delta _ s z _ i }$ の分布として、前回紹介した相乗対数正規モデルを用います。もちろん、このモデルが適用できない例もありますが、相乗対数正規モデルは良い近似になることが多いです。
　相乗対数正規モデルの非ガウス性は、非ガウスパラメタ $\lambda^ 2$ を以下の式を使って推定することで特徴づけられます。

$\displaystyle{ \begin{aligned} \hat{\lambda}^2_q &=\frac{2}{q(q-2)} \left[ \ln\!\left( \frac{\sqrt{\pi}\,\langle |\Delta _ s z _ i|^q\rangle}{2^{q/2}} \right)- \ln \Gamma \!\left(\frac{q+1}{2}\right) \right] & (q \neq 0, 2) \\ \hat{\lambda}_0^2 &=-\langle\ln | \Delta _ s z _ i| \rangle-\frac{\gamma+\ln 2}{2} &(q=0) \\ \hat{\lambda}_2^2 &=-\left\langle \Delta _ s z _ i^2 \ln \right|\Delta _ s z _ i| \rangle-\frac{\gamma+\ln 2}{2}-1 &(q=2) \end{aligned} }$

2. 独立同分布の非ガウス過程と間欠性過程の違い

　ここでは、数値的に生成した時系列の例を用いて、独立同分布（IID）の非ガウス過程と、間欠性をもつ過程の違いを具体的に見ていきます。以下のRスクリプトを使って生成される２つのサンプル時系列は、いずれも「非ガウス性」をもつ時系列ですが、その時間構造は本質的に異なります。興味があれば、以下のRスクリプトを実行してみてください。

############################################################
# IID 相乗対数正規過程 vs 間欠性（ランダムカスケード）過程
#  - 同じ長さ N (= 2^m) と同じ非ガウス性パラメタ λ^2 で生成
#  - 平均0・分散1に標準化して比較
#  - 分布はカーネル密度推定（density()）で推定し、同一グリッド上で重ね描き
#  - 分布図は裾を見やすくするため y 軸を対数表示
############################################################

############################################################
# 1) 生成関数
############################################################

#-----------------------------------------------------------
# IID 相乗対数正規過程： X = W * exp(Y)
#   W ~ N(0,1)
#   Y ~ N(-λ^2, λ^2)  （E[exp(Y)] = 1 となるよう平均を -λ^2 に設定）
#   ※分散は時間的に一様（間欠性なし）
#-----------------------------------------------------------
gen_IID_lognormal <- function(N, lambda2, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  w <- rnorm(N)
  y <- rnorm(N, mean = -lambda2, sd = sqrt(lambda2))
  x <- w * exp(y)
  return(x)
}

#-----------------------------------------------------------
# 間欠性（ランダムカスケード）相乗対数正規過程
#   - 重み（lognormal）を2分岐で m 段生成 → length = n0 * 2^m
#   - 最後に白色ガウス雑音を掛けて符号付き系列にする
#   - λ^2 を m 段に等分して総分散を揃える
#-----------------------------------------------------------
gen_intermittent_lognormal <- function(lambda2, m, n0 = 1L, seed = NULL) {
  stopifnot(is.numeric(lambda2), length(lambda2) == 1L, lambda2 >= 0)
  stopifnot(is.numeric(m),       length(m)       == 1L, m >= 1, m == as.integer(m))
  stopifnot(is.numeric(n0),      length(n0)      == 1L, n0 >= 1, n0 == as.integer(n0))

  if (!is.null(seed)) set.seed(seed)

  # 各段に割り当てる分散（合計が lambda2）
  lambda2_step <- lambda2 / m
  lambda_step  <- sqrt(lambda2_step)

  # 初期重み（すべて 1）
  weights <- rep.int(1, n0)

  # m 段の2分岐カスケード
  for (k in seq_len(m)) {
    n <- length(weights)

    # 2n 個の正規乱数 → 標準化 → 分散調整 → 平均シフト
    y <- rnorm(2L * n)
    y <- (y - mean(y)) / sd(y)
    y <- y * lambda_step - lambda2_step  # E[exp(Y)] = 1 を保つためのシフト

    # 2分岐（奇数・偶数に格納）
    new_weights <- numeric(2L * n)
    new_weights[seq.int(1L, 2L * n - 1L, by = 2L)] <- weights * exp(y[1L:n])
    new_weights[seq.int(2L, 2L * n,     by = 2L)] <- weights * exp(y[(n + 1L):(2L * n)])

    weights <- new_weights
  }

  # 最後に白色ガウス雑音を掛けて符号付き系列へ
  w <- rnorm(length(weights))
  x <- w * weights
  return(x)
}

############################################################
# 2) パラメタ設定（同じ N・同じ λ^2）
############################################################
lambda2 <- 0.6

# カスケード長を厳密一致させるため N = 2^m
m <- 12L
N <- 2^m

seed_iid <- 1
seed_int <- 1

# 図に表示する点数（先頭部分）
n_show <- 4000

############################################################
# 3) 時系列生成（標準化して比較）
############################################################
x_iid <- gen_IID_lognormal(N = N, lambda2 = lambda2, seed = seed_iid)
x_int <- gen_intermittent_lognormal(lambda2 = lambda2, m = m, n0 = 1L, seed = seed_int)

# 平均0・分散1に標準化（分布比較を見やすくする）
x_iid_z <- as.numeric(scale(x_iid))
x_int_z <- as.numeric(scale(x_int))

############################################################
# 4) カーネル密度推定（同じグリッドで評価）
############################################################

# 表示範囲：外れ値に引っ張られないように分位点で決める
q <- 0.999
x_lim <- range(
  quantile(x_iid_z, probs = c(1 - q, q)),
  quantile(x_int_z, probs = c(1 - q, q))
)

# 共通グリッド（同じ x 上で比較）
grid_x <- seq(x_lim[1], x_lim[2], length.out = 1000)

# density() を計算（内部グリッドは別でもOK）
dens_iid <- density(x_iid_z, n = 2048)
dens_int <- density(x_int_z, n = 2048)

# 共通グリッドへ補間して揃える
dens_iid_y <- approx(dens_iid$x, dens_iid$y, xout = grid_x, rule = 2)$y
dens_int_y <- approx(dens_int$x, dens_int$y, xout = grid_x, rule = 2)$y

############################################################
# 5) 作図（時系列2枚 + 密度1枚）
############################################################
par(mfrow = c(3, 1),
    mar = c(5, 5, 2, 1),
    las = 1, xaxs = "i",
    cex.lab = 1.5, cex.axis = 1.2)

#-----------------------------------------------------------
# (A) IID の時系列（先頭部分）
#-----------------------------------------------------------
plot(x_iid_z[1:n_show], type = "l",
     xlab = "i", ylab = "Standardized x[i]",
     col = 2,
     main = sprintf("IID 相乗対数正規過程（N=%d, lambda^2=%.2f）", N, lambda2))

#-----------------------------------------------------------
# (B) 間欠性過程の時系列（先頭部分）
#-----------------------------------------------------------
plot(x_int_z[1:n_show], type = "l",
     xlab = "i", ylab = "Standardized x[i]",
     col = 4,
     main = sprintf("間欠性（ランダムカスケード）相乗対数正規過程（N=%d, lambda^2=%.2f）", N, lambda2))

#-----------------------------------------------------------
# (C) 分布：カーネル密度推定（重ね描き、y軸は対数）
#-----------------------------------------------------------
ylim_y <- c(1e-4, 1)

plot(grid_x, dens_iid_y, type = "l",
     xlab = "Standardized x", ylab = "Density",
     col = 2, lwd = 2, yaxt="n",
     log = "y", ylim = ylim_y,
     main = "カーネル密度推定による分布比較（対数軸）")

# 主要目盛り（10^k）のみ表示（読みやすさ優先）
axis(2, at = 10^(-4:0), labels = parse(text = paste0("10^", -4:0)),
     las = 1, cex.axis = 1.2, tck = -0.02)

lines(grid_x, dens_int_y, lty = 2, col = 4, lwd = 2)

legend("topright",
       legend = c("IID", "Intermittent (cascade)"),
       lty = c(1, 2), lwd = 2,
       col = c(2, 4),
       bty = "n")

Rスクリプト（上記）の実行例。上：独立同分布（IID）の非ガウス過程、中：間欠性をもつ過程、下：両者の確率密度関数の比較。

　上のRスクリプトでは、２つの過程の条件を意図的に揃えています。

時系列の長さ N は同じ
非ガウス性の強さを表すパラメタ $\lambda^ 2$ も同じ
平均 0、分散 1 も同じ

つまり、同じ非ガウス性をもつにもかかわらず、時間構造だけが異なる状況を作っています。

2-1. IID 相乗対数正規過程とは何か

　最初に生成しているのは、IID な相乗対数正規過程です。この過程では、各時点の値は

ガウス雑音 $W$ （独立同分布過程）
対数正規的な振幅 $e^ Y$ （独立同分布過程）

の積として与えられます。確率分布は裾の厚い非ガウス分布になりますが、異なる時刻に出現する値は互いに独立です。

2-2. 間欠性（ランダムカスケード）過程とは何か

　間欠性を示す過程は、乱流研究で用いられてきたランダムカスケードモデルを参考にして生成したものです。このモデルでは、局所的な分散が階層的に分岐する形で相乗的に変調されます。具体的には、あるスケールで与えられた揺らぎ（標準偏差）の大きさが、より小さなスケールへと分割されながら確率的に（対数正規分布に従い）変化します。この操作によって、

大きな振幅をもつ区間が局所的に集中し
その周囲には比較的穏やかな区間が広がる

という特徴的な構造が自然に生じます（詳しい説明は参考文献を読んでください）。
　ここでは、非ガウス性の強さを表すパラメタ $\lambda^ 2$ を上のIID 過程と同じに設定しています。そのため、全体として見た非ガウス性の強さは同程度ですが、揺らぎの分散は時間的に強く不均一となり、明確な間欠性が現れます。

2-3. 時系列の違いと分布の一致

　上の図の上段・中段では、両過程の時系列を同じスケールで表示しています。

IID 過程では、大きな変動は単発的に散らばって現れます
間欠性過程では、大きな変動がバースト（塊）として現れます

この違いこそが、「非ガウスだが間欠性はない」過程と、「非ガウスかつ間欠的」な過程の決定的な違いです。
　一方、下段に示したカーネル密度推定による分布を見ると、両者は驚くほどよく似ています。特に対数軸で表示すると、裾の厚さもほぼ重なって見えることが分かります。これは重要なポイントです。「分布だけを見ても、間欠性の有無はほとんど区別できない」という事実を、この図は示しています。非ガウス性は「どれくらい大きな値が出やすいか」を表しますが、それがいつ・どのように現れるかまでは教えてくれません。

2-4. なぜ多重スケール確率分布解析が必要なのか

　このように、

単一スケールの確率分布はほぼ同じ
しかし時間構造はまったく異なる

という状況では、平均、分散、観測時系列の確率分布だけでは不十分です。もっといえば、上の例では両過程のパワースペクトルも同じになりますので、パワースペクトルも役に立ちません。そこで必要になるのが、粗視化スケールを変えながら確率分布の変化を追うという視点です。下の図が粗視化スケールを変化させた結果を描いたものです。

上図に見られるように、IID 過程では分布が $s = 32$ 程度で、急速に正規分布に収束するのに対し、間欠性過程では裾が厚く中心部が尖った非ガウス性が長いスケールまで残っています。目で見ると違いが良く分からないという人もいるかもしれませんが、赤と青の線で描かれた分布形状に何となく違うことはわかると思います。この違いを定量化するのが次の課題です。既に説明したように相乗対数正規モデルをフィットするのが一つの方法です。これについては、次の記事で解説することにします。
　粗視化時系列を生成し、分布を見積るRスクリプトを以下に掲載しておきます。

############################################################
# IID 相乗対数正規過程 と
# 間欠性（ランダムカスケード）相乗対数正規過程の比較
#
# 本スクリプトで行うこと：
#   1) 同じ長さ N (= 2^m)、同じ非ガウス性パラメタ λ^2 をもつ
#      2 種類の時系列を生成する
#   2) 粗視化スケール s を与え、スライディング窓で部分和を計算する
#        Δ_s z[i] = Σ_{j=0}^{s-1} x[i+j]   (i = 1, …, N-s+1)
#   3) 粗視化系列を平均 0・分散 1 に標準化する
#   4) 以下を作図する
#      - 粗視化後の時系列（IID / 間欠性）
#      - カーネル密度推定による分布（重ね描き）
#        ※ 裾を見やすくするため y 軸は対数表示
#        ※ 標準正規分布 N(0,1) を基準として重ねる
############################################################

############################################################
# 0) パラメタ設定
############################################################
lambda2  <- 0.6      # 非ガウス性パラメタ（対数振幅の分散 λ^2）
m        <- 15L      # カスケード段数（N = 2^m）
N        <- 2^m
s        <- 64L      # 粗視化スケール（窓幅）
seed_iid <- 3        # IID 用乱数シード
seed_int <- 3        # 間欠性過程用乱数シード
n_show   <- 10000    # 時系列プロットで表示する点数
q_clip   <- 0.999    # 分布表示範囲を決める分位点

############################################################
# 1) 時系列生成関数
############################################################

#-----------------------------------------------------------
# IID 相乗対数正規過程
#
#   X = W * exp(Y)
#     W ~ N(0,1)
#     Y ~ N(-λ^2, λ^2)
#
#   E[exp(Y)] = 1 となるように平均を -λ^2 に設定
#   → 分散は時間的に一様（間欠性なし）
#-----------------------------------------------------------
gen_IID_lognormal <- function(N, lambda2, seed = NULL) {
  stopifnot(is.numeric(N), length(N) == 1L, N >= 1, N == as.integer(N))
  stopifnot(is.numeric(lambda2), length(lambda2) == 1L, lambda2 >= 0)

  if (!is.null(seed)) set.seed(seed)

  w <- rnorm(N)
  y <- rnorm(N, mean = -lambda2, sd = sqrt(lambda2))
  x <- w * exp(y)
  x
}

#-----------------------------------------------------------
# 間欠性（ランダムカスケード）相乗対数正規過程
#
#   ・対数正規重みを 2 分岐で m 段生成
#   ・全体の対数振幅分散は λ^2（各段に λ^2/m を割り当て）
#   ・最後に白色ガウス雑音を掛けて符号付き系列にする
#-----------------------------------------------------------
gen_intermittent_lognormal <- function(lambda2, m, n0 = 1L, seed = NULL) {
  stopifnot(is.numeric(lambda2), length(lambda2) == 1L, lambda2 >= 0)
  stopifnot(is.numeric(m),       length(m)       == 1L, m >= 1, m == as.integer(m))
  stopifnot(is.numeric(n0),      length(n0)      == 1L, n0 >= 1, n0 == as.integer(n0))

  if (!is.null(seed)) set.seed(seed)

  # 各段に割り当てる分散
  lambda2_step <- lambda2 / m
  lambda_step  <- sqrt(lambda2_step)

  # 初期重み（すべて 1）
  weights <- rep.int(1, n0)

  # m 段のランダムカスケード
  for (k in seq_len(m)) {
    n <- length(weights)

    # 正規乱数を生成 → 標準化（数値安定性のため）
    y <- rnorm(2L * n)
    y <- (y - mean(y)) / sd(y)

    # 分散調整 + 平均シフト
    #   各段で E[exp(Y)] = 1 を保つ
    y <- y * lambda_step - lambda2_step

    # 2 分岐
    new_weights <- numeric(2L * n)
    new_weights[seq.int(1L, 2L * n - 1L, by = 2L)] <- weights * exp(y[1L:n])
    new_weights[seq.int(2L, 2L * n,     by = 2L)] <- weights * exp(y[(n + 1L):(2L * n)])
    weights <- new_weights
  }

  # 白色ガウス雑音を掛けて符号付き系列にする
  w <- rnorm(length(weights))
  x <- w * weights
  x
}

############################################################
# 2) 粗視化（スライディング窓による部分和）
############################################################
# Δ_s z[i] = Σ_{j=0}^{s-1} x[i+j]
# filter() を使って高速に計算
############################################################
coarse_grain_sum <- function(x, s) {
  stopifnot(is.numeric(s), length(s) == 1L, s >= 1, s == as.integer(s))
  if (s == 1L) return(as.numeric(x))

  # 右寄せの移動和（先頭 s-1 点は NA）
  cg <- stats::filter(x, rep.int(1, s), sides = 1, method = "convolution")
  cg <- as.numeric(cg)

  # NA を除去して長さ N-s+1 の系列にする
  cg[seq.int(s, length(cg))]
}

############################################################
# 3) 時系列生成 → 粗視化 → 標準化
############################################################
x_iid <- gen_IID_lognormal(N = N, lambda2 = lambda2, seed = seed_iid)
x_int <- gen_intermittent_lognormal(lambda2 = lambda2, m = m, n0 = 1L, seed = seed_int)

# 粗視化
dz_iid <- coarse_grain_sum(x_iid, s)
dz_int <- coarse_grain_sum(x_int, s)

# 平均 0・分散 1 に標準化
dz_iid_z <- as.numeric(scale(dz_iid))
dz_int_z <- as.numeric(scale(dz_int))

N_cg <- length(dz_iid_z)  # = N - s + 1

############################################################
# 4) カーネル密度推定（同一グリッドで比較）
############################################################
x_lim <- range(
  quantile(dz_iid_z, probs = c(1 - q_clip, q_clip)),
  quantile(dz_int_z, probs = c(1 - q_clip, q_clip))
)

x_lim_max <- max(abs(x_lim))
grid_x <- seq(-x_lim_max, x_lim_max, length.out = 1000)

dens_iid <- density(dz_iid_z, n = 2048)
dens_int <- density(dz_int_z, n = 2048)

dens_iid_y <- approx(dens_iid$x, dens_iid$y, xout = grid_x, rule = 2)$y
dens_int_y <- approx(dens_int$x, dens_int$y, xout = grid_x, rule = 2)$y

############################################################
# 5) 作図
############################################################
par(mfrow = c(3, 1),
    mar = c(5, 5, 2, 1),
    las = 1, xaxs = "i",
    cex.lab = 1.5, cex.axis = 1.2)

# (A) IID：粗視化時系列
plot(dz_iid_z[1:min(n_show, N_cg)], type = "l",
     xlab = "i", ylab = "Standardized Δ_s x[i]",
     col = 2,
     main = sprintf("IID（粗視化 s=%d, N=%d → %d, λ^2=%.2f）", s, N, N_cg, lambda2))

# (B) 間欠性：粗視化時系列
plot(dz_int_z[1:min(n_show, N_cg)], type = "l",
     xlab = "i", ylab = "Standardized Δ_s x[i]",
     col = 4,
     main = sprintf("間欠性カスケード（粗視化 s=%d, N=%d → %d, λ^2=%.2f）", s, N, N_cg, lambda2))

# (C) 分布：カーネル密度推定（対数軸）
ylim_y <- c(1e-4, 1)

plot(grid_x, dens_iid_y, type = "l",
     xlab = "Standardized Δ_s x", ylab = "Density",
     col = 2, lwd = 2, yaxt = "n",
     log = "y", ylim = ylim_y,
     main = sprintf("粗視化系列の分布（カーネル密度推定, s=%d）", s))

# 標準正規分布（基準）
curve(dnorm(x), add = TRUE, col = gray(0.2), lty = 2, lwd = 2)

axis(2, at = 10^(-4:0),
     labels = parse(text = paste0("10^", -4:0)),
     las = 1, cex.axis = 1.2, tck = -0.02)

lines(grid_x, dens_int_y, col = 4, lwd = 2)

legend("topright",
       legend = c("IID", "Intermittent (cascade)", "Gaussian (N(0,1))"),
       col    = c(2, 4, gray(0.2)),
       lty    = c(1, 1, 2),
       lwd    = c(2, 2, 2),
       bty    = "n")

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。